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SYNOPSIS 


In many branches of Science and Engineering applications / 
one often encounters singular boundary value problems in ordinary 
differential equations, which although mathematically well 
conditioned, are virtually difficult to solve with analytical 
procedures* For solving such problems, numerical methods have 
received a great deal of attention in the last few years* 
Basically, there are two approaches to solve these problems/ 
one is to apply directly the numerical methods to singular 
problems ignoring the singularity and second is to give some 
analytical treatment to singularity and then solve the resulting 
regular boundary value problem* In the former approach, one 
would generally need a very fine mesh especially in the vicinity 
of the singularity thereby increasing the computational complexity 
of the problem* One may also sometimes need to consider higher 
order methods to achieve a reasonably good accuracy* In the 
latter case, one expects to produce good results in view of the 
analytical treatment to remove the singularity in the problem* 
However, the ingenuity of the approach used to subtract out 
the singularity and then to solve the resulting regular problem, 
would determine the efficiency of the method* 

In the present thesis, we have derived some new and 
efficient computational algorithms for solving singular two 
point boundary value problems over finite and infinite intervals* 

Chapter 1 presents introduction, types of singularities, 
motivation with physical applications, a brief survey of the work 
done by earlier researchers and summary of the work done in the 
present thesis* 



In -the second chapter, a particular class of singular 
boundary value problems have been considered for ordinary 
differential equations on a finite interval, which has a 
singularity of the first kind* These types of problems have 
been discussed earlier by researchers,,and arise when reducing 
partial differential equations by physical symmetry* The 
singularity is subtracted out by using the series expansion 
and a new boundary condition is derived in the vicinity of 
the singularity# A discrete invariant imbedding method .has 
been developed to solve the resulting regular boundary value 
problem efficiently# The stability of the method is discussed# 

The proposed method is conceptually simple and easy to implement. 
Some test examples have been solved and the numerical results 
presented# A comparison of numerical results has also been 
made with other methods# 

Chapter 3 deals with the continuous version of Invariant 
Imbedding method to solve the singular boundary value problems 
considered in Chapter 2# Again the expansion procedure has been 
employed to derive the new boundary condition near to singularity* 
Next a continuous invariant imbedding has been described to 
reduce the regular boundary value problem into a set of initial 
value problems* The Runge Kutta Pehlberg variable step size 
technique has been employed to integrate the stable initial 
value problems efficiently* Some model problems have been 
solved to illustrate the efficiency of the method* The continuous 



iii 


Invariant Imbedding combined with variable step size initial 
value software have been found to be more accurate than the 
discrete version of Invariant Imbedding but at a slightly more 
cost of computational time* 

The treatment given to the singularities to the problems 
considered in earlier chapters is based on using extended power 
series and thereby derive a new boundary condition near to 
singularity* However, in cases where there is slow convergence 
of the series expansion, an economized expansion may be used* 

In Chapter 4, we discuss Chebyshev economization near the 
singular point to overcome slow convergence of the series* 

Using the shifted Chebyshev polynomials and the Lanczos Tau 
method a new boundary condition near the singularity is 
obtained* The resulting regular boundary value problem is 
then treated with continuous Invariant Imbedding described in 
Chapter 3* Numerical results of some test examples are also 
included* 

The attempts in earlier chapters for the removal 
of singularity are based on using expansion procedures in the 
neighbourhood of the singularity* However for some problems, 
it may be difficult or even not possible to obtain the series 
expansion near the singularity* In such cases, one may have to 
use a direct approach to solve singular boundary value problems* 

A modified fourth order finite difference method to solve a 
certain class of linear singular BYP^s is presented in Chapter 5* 



The original differential equation is modified at the singular 
point* The main feature of the modified difference scheme is 
that it leads to the tridiagonal system of equations which has 
been solved by discrete Invariant Imbedding method* Some model 
problems have been solved and the numerical results included* 

The numerical solution of nonlinear singular boundary 
value problems has been discussed in Chapter 6* A quasi- 
linearization technique has been used to reduce the nonlinear 
problem to a sequence of linear problems* A fourth order finite 
difference method described in Chapter 5 has been employed to 
solve these problems effectively* Some model problems have been 
solved to demonstrate the efficiency of the method* 

In Chapter 7, we have considered the infinite interval 
boundary value problems in ordinary differential equations • 

In Case I, We derive the asymptotic boundary condition to 
represent the infinity condition at the far end* This boundary 
condition has been derived such that it expresses the asymptotic 
behaviour of the solution as well, and converge to the actual 
solution of the * infinite 7 problem as the length of finite 
interval tends to infinity* The resulting boundary value problem 
is treated by imbedding methods * The stability of the problem 
is analysed* In Case II, we have briefly included the 
transformation technique as an alternative method to solve 
these problems* For the transformed problems on a finite 
interval with a regular singularity at one end point of the 
interval, the methods discussed in earlier chapters can be 
applied directly* 



V 


Finally in Chapter 8, the application of cubic spline 
method to infinite interval problems is discussed* By reducing 
the infinite interval to a finite interval which is large and 
imposing the asymptotic boundary conditions at the far end, the 
two point boundary value problem on a finite interval is 
treated by using cubic spline approximation* The tridiagonal 
system resulting from the spline approximation is efficiently 
solved by the method of imbedding* The stability of the method 
is also discussed and the theory is illustrated by test examples* 

In a nut- shell, the numerical methods presented in this 
thesis for solving singular boundary value problems over finite 
and infinite intervals have been shown to be efficient over other 
methods* Although problems with singularity at one end point 
have been discussed, yet one may apply these methods to problems 
with interior regular singularities or two regular singular end 
points* The invariant imbedding methods is conceptually simple, 
easy to implement, and has the properties of solving boundary 
value problems subject to different boundary conditions without 
extra computational effort* 

Model problems have been solved and numerical results 
presented in respective chapters * The numerical results presented, 
in this thesis have been computed on DEC— 1090 system in single 
precision arithmetic* 



INTRODUCTION 


In the recent past/ methods for numerically solving two 
point boundary value problems for differential equation have 
been extensively studied by several researchers* A number 
of methods like Shooting Methods [l7 , 55 , 78 ] , Finite 
difference methods [27 , 55 , 57 / 63 / 67 , 81]/ Collocation 
Methods [l* 18/ 82] and Invariant imbedding methods [6 * 11 / 

13 / 15 / 19 / 58 * 69-71/ 11 j 19 t 85-90] have been reported 
in the literature for solving these problems* However little 
effort seems to have been devoted to the solution of Singular 
two point boundary value problems in ordinary differential 
equations* Difficulty arises in solving Singular boundary 
value problem due to the fact that the coefficient functions 
in the differential equation are not analytic* When the 
analytical methods fail to produce the solution of Singular 
boundary value problem (SBVP)/ one then should try to ascertain 
the approximate nature of the solution* Therefore/ it is very 
useful to develop numerical methods to obtain approximate 
solution to hard problems like SBVP* Such attempts may prove 
to be useful in many practical situations where/ analytical 
solutions are difficult to obtain* The purpose of this thesis 
is to devise and develop efficient computational methods for 
solving singular boundary value problems over finite and 
infinite Intervals* It would be pertinent here to first give 
the classification of singularities and their occurrence in 
practical applications* 
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1*1 Types of Singularities 

We begin the process of local analysis by classifying 
the singularities. Singularities can be of two kinds : 

( i) Singularity of the first kind which is also known as 
Regular Singularity, and ( ii) Singularity of the second kind, 
which is sometimes refered to as essential singularity* We 
classify a point x Q , which may be an Ordinary point, a regular 
singular point or an irregular singular point of a differential 
equation as follows : 

Suppose we consider a second order linear differential 
equation of the form, 

f 0 (x)y"(x) +f^(x)y / (x)+f 2 (x)y(x) = f 3 (x), x e[ x Q ,b] (l#l) 


where f 0 (x), f^(x) and f 2 (x) are analytic in the neighbourhood 
of some point x = x^say), then x = x^ is said to be an 
Ordinary point in the sense that the solution of (1*1) can be 


represented by a Taylor series in powers of (x—Xq) * If for 
some point x = Xq, f 0 (x^) = 0 and f^Cx^ 0, then Xq is called 
a Singular point of the equation ( 1 .1 ) . In such a case 
rewriting Eqn* (l*l) in the form 


y"( x) +F^ ( x) y * ( x) +P 2 1 x) = F 3 (x)' ( 1 *2) 


where F^( x) 
F^(x) and F 


f , ( x) 

= , i = 1,2, we see that the coefficients 

2 (x) fail to be analytic at x = x Q * 
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The point x = is said to be a Regular ^Singular- joint 
of (1*2) if (»-x 0 )F^(x) and ( xt-Xq) ** 2 ^ are analytic at- x =» 3^* 

A solution of (1*2) may be analytic at a regular singular point* 
If it is not analytic its singularities must be either a pole or 
an algebraic or logarithmic branch point. It is pointed that 
there is always one solution of the form 

y(x) = (x-x^) a A(x) 

where a is a number called the indicial exponent and A(x) is a 
function which is analytic at and which has a Taylor series 
whose radius of convergence is atleast as large as the distance 
to the nearest singularity of the coefficient functions* 

The point x = x^ (Xq £ *>) is called an Irregular S^ingular 
point of (1*2) if it is neither an ordinary point nor a regular 
singular point* There is no comprehensive theory of irregular 
singular points, we can say that at an irregular point, all 
functions exhibit as essential singularity, often in combination 
with a pole or an algebraic or logarithmic branch point* 

Often it is of interest to investigate solutions of the 
form (1*2) for large values of x, say x « «•* A simple way of 
doing this is to make use of the substitution x = jg in (1*2), 
and study the solutions of the resulting equation near t =* 0# 

Then, for example, the results on analytic equations and equations 
with a regular singular point t = 0 can be applied* Using the 
substitution x «* ^ Eqn* (1*2) reduces to the form 
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t 2 y« + (2 


F.,(t) P 2 (t) F-(t) 

t ■*-) ty* + * *y y *= " 2 " 

t t 


(1.3) 


where F^(t) ' = F^(^), P^Ct:) = F 2 (:£) F 3 (t) = F 3^^* We sa Y 

that infinity is a regular singular point for (1*2) if the 
Eqn* (1*3) has the origin t = 0 as a regular singular point. 

If t = 0 is not a regular singular point, then infinity is said 


to be an irregular singular point* 

The simplest example of an equation with a regular 
singular point at infinity is the Euler's equation* viz* 


2 

a y** + uxy 7 + by = 0 

where a,b are constants* This equation has the origin and 
infinity as regular singular points* 

1*2 Motivation and Physical Problems : 


Singular boundary value problems for ordinary differential 
equations arise very frequently in several areas of science 
and Engineering* For example, consider the following boundary 
value problem 


d 2 T , P 
dx x 


dT 


+ f(T) 



dT(0) 

dx 


0, T(l) = 1 


(1.4) 


which results from an analysis of heat conduction through a 
solid with heat generation* The function f(T) represents the 
heat generation within the solids this, in general, is a 
function of the local temperature T* The constant p is equal 
to 0, l or 2 depending on whether the solid is a plate, a 
cylinder or a sphere* In a simple case when p = 0, the solid 
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is a flat plate with constant heat generation and f(T) a q/K, 
where K is the heat conductivity and q is a heat generation per 
unit volume* Eqn. (1*4) then becomes 



n dT(O) 
' dx 


0, T(l) = 1 


which is a linear differential equation with constant coefficients* 
Thus, we arrive at a two point boundary value problem ■which 
has been solved efficiently by many authors* 

Consider next the case in which the solid is a cylindrical 

rod and the heat generation is linearly proportional to the 

2 

temperature T. Then f(T) = (3 T with p = 1, and Eqn* (1*4) 
becomes 


d 2 T 1 

: '2 + - 

dx 


_ ^T 
x dx 


+ 0 T = 0 


dT(0) 

"dx ’ 


d(l) = 1 (1.5) 


which is a linear differential equation with variable coefficients* 
For these types of problems the solution becomes complicated due 
to the singularity arising in the differential equation. One 
then needs to take extra care to solve these problems* Similarly 
when p = 2 and the heat generation is written as f(T) = ae , 
where a is a constant, Eqn* (1 *4) becomes 


d 2 T ^ 2 
dx 2 + * 


dT 

d^ 


T 

+ ae = o. 


dT(0\ 

dx 


= 0,. T(l) = 1- 


( 1 . 6 ) 


This problem becomes a nonlinear boundary value problem due to 
the nonlinearity of the heat generation term. For such types 
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of problems, the numerical methods seem, to be the only choice 
for their solution* 

To quote a few more examples of practical importance 
in different areas, we have : 

i) In Astronomy : The equilibrium of isothermal gas spheres 
can be described by 


y"(x) + ^ y* (x) + (y(x)) 5 = o 
y'(0) = 0 / y(l) = ^3 

This problem is of the form (l*4)« 


(1*7) 


ii) In Chemistry : The formulation of heat and mass transfer 
within porous catalyst particle leads to [48] 


_ ~ rp ( l-y( x) ) 

y"(x) +“ y r (x) a* y(x) exp [• — — -] 

l+p(l~y(x)) 

y f ( o) = O? y(l) = 1 


( 1 . 8 ) 


where 

y : dimensionless concentration, 

a : Structure of Catalyst particle, a = 2 spherical, 
: Chemical constants* 


iii) Thomas— Fermi Model : The Thomas— Fermi model in atomic physics 
describes the charge concentration y( x) of electrons in 
an ion : 
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y"(x) = x 1 (y(x)) 1/2 

(1-9) 

y(0) = i ; y(b) =o* 

iv) The Ginzburg-Landau Equations : A problem in superconductivity 
to find cylindrically symmetric solutions of the Ginzburg— 
Landau equations with radius r, fluxoid quantum number n = 1/ 
and G in zbu rg-Lan dau parameter K = 1 leads to 


L [f(x) ] = K 2 f(x) (f 2 (x) - 1 + a 2 (x) - ^ a(x)) 
L [a(x) ] = f 2 ( x) (a(x) - 


( 1 - 10 ) 


a(0) * 0 ; a(r) = R 

0 < X < r 


f(o) =0 ; f'(r) = 0 

where f(x) : order parameter, 
a(x) : vector potential, 
R : free parameter. 


r 

and L [ • ] 


« 10 

_ ah.) .1 a(.) 

ax 2 x dx 




v) Thin Shallow Spherical Shell : The elastic stability o'f 
■ thin shallow spherical shells subject to uniform pressure 
is described by [ 67,9 1] 


f"(x) = - ju 2 g(x) - 2y + f(x) g(x) - ^ f*(x) 
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g /y (x) = w 2 f(x) - f 2 (x) ~ ^ g*(x) 

f'CoJ = g'(o) = 0/ Cl.li) 

f(l) = O t g y (l) 4- ( 1-y) g( 1 5 =0 clamped edge condition* 
where 

x : normalized polar angle 

f(x) s normalized angular deflection 

g(x) : normalized stress 

V : load parameter 

2 

1± : geometry of the shell 

v = (Poisson's ratio). 

A few notable examples of boundary value problems over 
infinite intervals are the Von Karman Swirling flows [61 J* a 
Combined forced and free convection flow over a horizontal 
plate [84], an eigen value problem for the Schrodinger 
equation [6 2} , Falker-Skan equation in boundary layer theory [67]* 
Unsteady flow of power-law fluids [67]* Flow through a porous 
medium [67] and several others. 

In the next section* we give the survey of literature 
on regular singular boundary value problems over finite intervals. 
In section 1*4* the work on infinite interval boundary value 
problems* done by several researchers is given* The work* done 
on boundary value problems with essential singularity is given 
in section 1*5. Finally, a short summary of the work done in the 
present thesis is included. 
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1*3 Singular boundary value problems with a Regular singularity : 

F or the last two decades, considerable number of 
methods have been discussed for solving Singular boundary value 
problems in ordinary differential equation with a singularity 
of first kind on a finite interval*. These methods have been 
developed either to treat singular boundary value problems 
directly or to first give local treatment to singularity and 
then solve the resulting regular problem* We begin with a 
method based on Finite difference analysis* A numerical method 
for generalized axially symmetric potential of a function 
u(x, y) described by 

L k [u] = AU + = 0 (1.12); 

satisfying certain boundary conditions was first studied by 
Parter [72] * However this particular type of partial differential 
equation can be reduced to the differential equation of the form 

2 

^ ^ + x Hx ^ = Jki < 1, q > 0 (l*13) 

dx' 

by separation of variables for the equation (l*12) in a rectangxe* 
Jamet [50] considered a linear ordinary differential equation 
of the second order ' 

Lu = ^ + f(x) — g(x)u = H(x) (1*14) 

dx^ 
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where 

f(x) e C(0,l] 
f(x) ■* ® as x •* 0 
g(x), H(x) e c[o.,l] 
g(x) > 0. 

According to rate of growth of f(x) near the origin, the problem 
to be posed may be a two point boundary value problem or one 
point boundary value problem* The author studied the application 
of three point finite difference approximation with a uniform 
mesh h and showed that the error is 0(h^ _cr ) when f(x) < — , for 
x small , a a constant with 0 < a < 1* The existance and 
uniqueness of the solution alongwith the convergence of certain 
finite difference schemes have been discussed* 

In another work, Ciarlet [24] has considered a nonlinear 
two point boundary value problem whose coefficients have a 
singularity of first kind* Based on Ritz— Galerkin' s piecewise 
linear approximation, the error in the uniform norm of 0(h ) 

has been given* This work is the extension of the authors work 
for a nonlinear two point boundary value problem [23] « 

Application of finite difference method to singular 
problems has also been discussed by Gustaffsson [47 ] for a scalar 
case* IXie to effect of singularity, if a difference approximation 
is applied on the whole interval, the convergence rate will be 
very poor* To improve upon the convergence, a series solution 
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can be represented on a subinterval near the singularity* to 
obtain a regular boundary value problem on the remaining interval * 
A centered rth order approximation is used to get a convergence 
rate 0(h ) on the whole interval* provided the series solutions 
are calculated with rth order accuracy* 

An approach via generalized projection method involving 
appropriate generalized spline functions which prove the converge 
faster than the usual finite difference scheme has been used by 
Natterer [68] for the numerical treatment of first order system 
of ordinary differential equation with a regular singularity* 

The author considers the approximation of the solution of 
singular boundary value problems such as 

y" + x~%y + Wy * f, 0 < x < 1* (1.15) 

where y is a vector in R n , U is a constant matrix and where the 
matrix W and the vector f depend continuously on x* In fact 
more general forms are considered where a similar singularity 
is allowed at the other end point. With the differential 
equation (l*15) are associated appropriate boundary conditions, 
the form of which depends upon the matrix U* After a careful 
study of the continuous problem, a projection method using 
generalized vector valued splines is described* The main result 
is that, with special choices of non-uniform meshes the 
convergence is the same as in nonsingular case up to a logarithmic 


f actor. 
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For a second order linear differential equation with a 
regular singular point/ Cohen and Jones [26] have applied an 
economized chebyshev expansion on the interval (0, l] instead 
of the series expansion* A second order finite difference 
method with deferred correction has been applied outside the 
range of economized expansion* 

The studies of more general linear and nonlinear 
singular boundary value problems have been made by Russel and 
Shampine [83]* The aim of the study is to compare the 
applicability of the Finite difference. Collocation and Patch 
bases procedures and the ease and effectiveness of their 
numerical comoutation. The collocation with piecewise 
polynomial functions in different sub intervals with a particular 
choice of collocation points for nonlinear two point boundary 
value problems has been described* A three point finite 
difference methods along the lines of Jamet [50] and Rose [81] 
has been used in their finite difference analysis* Patch bases 
method discussed earlier by Rose [8l] has been extended to 
Singular boundary value problems* The authors feel that in many 
circumstances one ought to use the traditional finite difference 
methods* The collocation results presume the mesh space h is 
sufficiently small, to even guarantee the approximate solution 
exist* They emphasized one might use the collocation if higher 
order procedures of uniform mesh or if one has complicated 
boundary conditions* 
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F *R* de Hoog and Weiss [ 29] have studied the solution for 
a first order system with a singularity of first kind. The 
application of some standard finite difference schemes (Box 
scheme and Euler's scheme) for the problem 

Y* ( t) - -jr y( t) a f ( t, y ) , 0 < t < 1 

with (1*16) 

b(y(0) , y(l) ) = 0 

and the linear eigen value problem 

y'Ct) - y(t)-A( t)y( t) = AG(t)y, 0 < t < 1 

with (1*17) 

b(y(0) / y(l) ) = 0 

with necessary boundary conditions, has been investigated* The 

author' s restrict their attention to the case when the solutions 

to these problems are continuous on [ 0, l] and differentiable 

on (0,1], The main aim of this work is to derive the most 

general boundary condition, which yields a Fredholm alternative 

for the problems* Under natural assumptions the finite difference 

schemes are shown to converge and have the usual convergence 

(up to possible logarithmic term for box scheme) provided 

y, f and b are sufficiently smooth* The asymptotic expansion 

of the error is then examined for a number of specific cases and 

2 

shown that the Keller's Box scheme has an 0(h ) expansion for 
many important problems. 
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In a recent -work* Reddien [73] studied a collocation 
method for the numerical solution of singular boundary value 
problems* These methods were developed by considering Certain 
projection on to a finite dimensional linear spaces of singular 
polynomial splines* The purpose of this note is to enhance the 
applicability of the methods by showing that the class of singular 
splines used in [74] possess convenient local support bases 
which are of considerable advantage in the actual compu tat ions* 

A detailed description of a shooting algorithm based 
on a Taylor series method is discussed by Rentrop [75] for the 
solution of a two point boundary value problem of the form 

y # = fCt^y), 

r N -■ N N /v 

y : [a^b] •* R f f : [a,b] x R -* R Cl *18) 

r(y(a)/y(b)) = 0; r : R N x R N -*■ R N 

The author extends this work to singular problem# The algorithm 
computes a Taylor series expansions in a small interval in order 
to kill the numerical singularity and a multiple shooting 
technique is used to solve the problem on the remaining interval* 
The reliability of the method is demonstrated by solving the 
G inzburg-Landau equations [76] arising in the theory of 
superconductivity* 

A numerical method for linear systems of first order 
equations with a regular singular point at one end point has been 
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given by Brabston and Keller [l6] * A class of problems 
considered for their analysis is, 

Ly(t) = y' ( t) - A(t)y(t) = f(t), 0 < t < 1 ; 

Ey(t) = lim [ B ft ( t)y( t)+B 1 y(-t)-b(t)l = 0 / (l»19) 

■ t*0 ° x 

A( t) = t~ 1 R + A q ( t) - 

where y(t), f(t), b(t) are n vectors while R,A Q ( t),B 0 ( t),B^ are 

nxn matrices ♦ The author's assumed A Q (t) Analytic on (0,6 Q ] 

and sufficiently smooth on (0, l] » Similarly B 0 (t),f(t) and 

b(t) are smooth on (0,l] tut may be singular at t * 0* The 

standard procedure of expanding about the singularity to get 

a nonsingular problem over a reduced interval is justified in 

some detail* The truncated regular two point problem is 

approximated by a stable rth order difference scheme (or 

centered Euler) as presented in [57] * The net is chosen to 

be uniform* The results computed are of 0(h ) accurate 

approximations* The author's also used Richardson extrapolation 
4\ 

to get 0(h) approximations# 

Jesperson [5l] studied the application of Ritz-Galerkin 
method for singular boundary value problems of the type arising 
when Poisson's equation -Au = f is encountered on a domain with 
cylindrical or spherical symmetry and is reduced to a one 
dimensional problem* The idea of the method is to derive a 
periori and norm estimates for the error. The difficulty 
is that these norms are not natural norms for the reduced problem# 
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With the aid of B— splines [18]* they proved some theoretical 
results and used these to derive the desired error estimates* 

In another work, de Hoog and Weiss [30] have analyzed 
the application of collocation methods to singular problems 
such as one considered by the same authors in [29] * They 
showed that collocation gives the same results for singular 
boundary value problems as for non-singular problems* 

Elder [38] showed how to apply the subscheme of 
invariant imbedding known as integration to blow up to confute 
the smallest eigenlength of a linear first order systems having 
a singularity of first kind* The purpose of this paper is to 
define an invariant imbedding for the singular TPBVP and then 
give the appropriate initial conditions at the origin for the 
Ricatti differential equation* They showed that the smallest 
eigenlength of the singular TPBVP is the right endpoint of the 
maximal interval of existence of the Ricatti initial- value 
problems* With this result it is reasonable to consider 
computing this smallest eigenlength by integrating the Ricatti 
equation to its singularity or to "blow up"* This method is 
referred to as " integration- to-blow up" in the papers by Boland 
and Nelson [15] and Nelson and Elder [70] • This procedure has 
also used in [l5,?o] with great success in the computation of 
eigenlengths of nonsingular TPBVP* 

Nelson, Sagong and Elder [7l] have considered a 
linear homogeneous first-order differential system with a 
singularity of the first kind at z = 0, 
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A B 

U r (z) a» (-^ + A^(z))u(z) + (-^ 4- B^(z))v(z}/ 

v'(z) = (-^ + C 1 (z))u(z) + (-— + D 1 (z))v(z)* 
Z 1 z 1 


( 1 * 20 ) 


where u and v are respectively in m- vector and ah n— vector of 
dependent variables/ and the matrices for i = 0/1 

are of the appropriate dimensions (e*g* is mxra) # The matrices 
W C o' and D Q are constant/ and A^ # B^ y C^, are analytic 

on (O/fi] for some <5 > 0 and piecewise continuous on [o,«>), with 
the boundary conditions of the form 


u( 0 +) and v( 0 +) exist (finite) 

and (l. 2 l) 

v(x) = 3 , x > 0 . 


They showed how the approach of Elder can be adapted to the 
solution of (1*20) - (l*2l) by means of Scott's version [85] of 
invariant imbedding* 

Recently Chawla [2l] discussed the construction of three 
point finite difference approximation and their convergence unaer 
appropriate conditions for the class of singular BVP's of the 
form 

(x <X y / )* = f(x/y)/ 

. . . I 

y(0) = A; 0 < a < 1 . (1*22) 


y(l) = B/ 
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The author established a certain integral identity on a general 
mesh from which the various methods are derived# For a non- 
uniform mesh over [0,1] / a method based on an evaluation of f 1 *- 

is obtained# For a uniform mesh/ two methods based on three 
evaluations of f are obtained* When a = 0, the method reduces 
to a classical second order method based on one evaluation of f* 
However 0(h ) convergence for all the defined methods is 
established# 

More recently Doedel [ 35/ 36] have constructed some high 
order finite difference methods for non-singular two point 
boundary value problems* Doedel and Reddien [37] have extended 
these methods to an important class of singular boundary value 
problems considered by Jamet [50] and other authors- A compact 
finite difference approximation the the differential equation 
have been considered and these have been shown to be equivalent 
to the collocation scheme# 

1*4 Boundary value problems on Infinite intervals : 

Another type of singular problems is the study of 
solutions on unbounded interval/ if the interval of interest is 
infinite or semi— infinite, where infinity is a regular singular 
point# For instance, the interval 0 < x < 00 sometimes called 
a semi— inf inite interval, since it does have one finite end 
point, a boundary condition would normally be imposed at x = 0, 
due to no boundary exists on the other end* These problems arise 
in many branches of science and Engineering* Problems of 
practical importance frequently occur in Fluid dynamics. 
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Aerodynamics/ Quantum Mechanics etcetera* 

Often in most cases/ the analytical solutions for these 
problems are not readily attainable and thus the problem is 
brought to the problem of finding efficient computational 
algorithms for obtaining numerical solution* Numerical methods 
for Infinite or semi-inf inite interval problems have been 
discussed by some authors* One of the ways of solving these 
problems is based on truncating the unbounded domain to a 
large but finite region. But one must make the decision as 
to how large the domain must be in order to properly represent 
an infinite or semi- inf inite domain in the computational sense. 

A method for solving linear two— point boundary value problem 
on an infinite interval 


L[y] = y" + p(x)y / + q(x)y = f(x) 
y(a) = a / y(oo) — o 


(1*23) 


is given by Pox [42/43/44]* Applying the second boundary 
condition at a finite point x = N, N arbitrary/ the solutions 
are computed by taking different values of N and observing the 
variations in the solutions till the desired accuracy is achieved* 


Robertson [80] has discussed a second order finite 
difference scheme for solving the Infinite interval problem 
given by (1*23) • The objective of his paper is to obtain 
numerical solution in a relatively small finite interval* The 
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method proposed is to examine y^( x) , the solution of Ly = f 
with y(a) = a, y(b^^) =0 under conditions which ensure 
y ( x) -* y(x) as N -*• <»♦ The results are shown to have O(h^) 
convergence# Several authors [ 49 , 67] adapted the procedure 
of choosing a large point to represent infinity before computing 
the numerical solutions# 

In another work, Alspaugh [5] studied the application 
of Invariant imbedding to the solution of boundary value problem 
on infinite interval* By using invariant imbedding, the author 
converts the boundary value problem to initial value problems 
and the resulting Ricatti's equations are integrated numerically 
until the desired accuracy is obtained# For certain ill- 
conditioned problems, the usual methods of solution of two point 
boundary value problems fail- It has been shown that the 
numerical stability of invariant imbedding formulation permits 
the easy solution of such problems* Several criteria for 
determining the appropriate length of integration are presented# 

Another way of solving the infinite interval problems is 
to reformulate the boundary value problem in a standard form# In 
such cases, one may use the coordinate transformations 
(Algebric or exponential transformation) to reduce the given 
problem to the problem over a finite interval# This procedure 
has been adopted earlier by several authors# Recently, Grosch 
and Orszag [46] investigated the utility of transformation 
techniques to solve boundary value problems in infinite regions# 
The utility of transformation has been shown through some 
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physical examples. However the great disadvantage of this 
approach is that the transformation technique usually produces 
a singular problem, which has an irregular singularity. 

Solutions that vanish rapidly or approach a constant at 
infinity are readily treated by mapping, but solutions that 
oscillate out to infinity are not so amenable to these techniques. 
The obvious advantage of these transformation is that no 
experimental choice of a large number, which represents the 
point at infinity is necessary. 

In a recent survey, Ascher and Russell [8] have shown 
how various non-standard problems can bo transformed to fit a 
simple format. Suppose a single equation expressing a boundary 
condition involves the value of the solution of the differential 
equation at more than one points it is shown how the problem 
can be formulated so that each boundary condition involves only 
one of the boundary points. Suppose a oroblem involves a 
boundary condition at an unknown point, it is shorn how to 
obtain an equivalent problem with known boundary points. Also 
reformulations to remove singularities and infinite intervals 
are discussed. 

Recently boundary value problems over semi-infinite 
intervals have been studied by Lentini and Keller [61-62], The 
authors deal with two point boundary value problems of the 
following forms 
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(a) 

y'( t) = A( t)y{ t) + f(t). 

t Q < t < «> ; with conditions 

(b) 

sup I l y( t) i I < co; 

^o 

(1-24) 

(c) 

B oY<t 0 ) + lim B oo y( t) = 3 
t-*oo 

( 3 = a constant vector) . / 

(a') 

= f(t ' u(t)) ' fc o ^ 

t < ■» * with conditions 

(b') 

lim u(t) = Uoo exists; 

t“^oo 

(1*25) 

( c' ) 

b 0 (u(t 0 ),uj = 0 , 



under suitable smoothness conditions for the matrix A(t) and the 
vector f(t/u(t)) and for appropriate constant matrices B o , 
and vector b Q - The technique of the reduction to a finite 
interval problem is adopted and its correct implementation is 
shovel# The basic idea employs that for linear problems the 
space S B of bounded solutions is finite dimensional and is 
easily determined- Then the proper 'boundary condition at 
infinity' is simply that the projection of the solution P n to 
the. complement of S B must vanish- Existence, uniqueness and 
behaviour of the solutions at infinity are discussed- The 
appropriate projection condition is simply applied to a finite 
point to obtain a finite domain- However a priori estimates 
for the location of finite boundary are not easy to obtain- 
In another work, the application of the technique has been 
shown by the same author for the Von Karman Swirling flows [61] , 
and a linear elasticity problem [60] - , 
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Recently a complete treatment, which give the asymptotic 
behaviour of the solution has been presented by some authors# 

F *R» De Hoog and Weiss [32] solved an infinite interval problem 
by restricting the problem to a large but Finite interval and 
imposing certain supplementary conditions at the far end# The 
success of this procedure depends upon the proper choice of 
these conditions* For a rather general class of problems the 
authors give a characterization of all possible supplementary 
boundary condition that work and examine the rate of convergence 
of the finite problem to that of original problem as the 
interval length of the finite problem tends to infinity and 
describe the boundary condition for which this rate is optimal# 

Markowich [64] formulated a adhoc method to solve 
boundary value problem on infinite interval of the form 


y' = t a f(t,y), 1 < t < »/ as 3N 


(1.26) 


y s C([i,co]) <==£> y e C [ljfoo)) and lim y(t) exists 

t-+oo 


where f : 


b(y(l) ) = 0 


where f : R n+ ^ -*■ R n and N Q is the set of nonnegative integers# 
The author deals with the numerical solution of boundary value 


problem by truncating the infinite interval to a finite but 
large one and to impose additional boundary condition at the for 
end* These boundary conditions should be posed in a way so that 
they express the asymptotic behaviour of the actual solution# 
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The resulting boundary value problem can be reduced the 
following form : 

X* T = t^tyXp), 1 < t < T, t >> 1, 

b( XrJ,( 1 ) ) =0 (1*27 ) 

sCxjCt^t) = 0 

In another paper, Markowich [65] studied the solution 
of the system 

y* = t a f(t,y), 1 <‘ t < «>, a > 0 (1*28.) 

with separated boundary conditions at 1 and <» where appropria te 
asymptotic boundary conditions are determined at °° to ensure a 
smooth solution* Asymptotic analysis of solution behaviour 
is used to determine the point T = T(e) such that the solution 
to a regular BVP on [l,T] is within 0(e) of the original 
solution* Certain symmetric (spline) collocation schemes are 
analyzed for solving this BVP* By equidistributing truncation 
errors, meshes with exponentially increasing subinterval 
lengths, the solutions with 0(e) accuracy are produced* The 
authors also show that the condition number of the collocation 
equations is asymptotically proportional to the number of 
mesh points employed when using this exponentially graded mesh* 

For K Gauss points collocation, this is only 0(e » 

Stability as e ■* 0 ( T(s ) -* 00 ) is shown for the linear case* 
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1*5 Essential Singularities : 

Recently the numerical solution of boundary value 
problems with an essential singularity has been investigated 
by Hoog and Weiss [31] * The boundary value problems of the 
type 


t^V = f(t,y), 0 < t < 1 , 

b(y(0),y(l)) = 0 


y e c[o,l] O c' ( 0 , 1 ] 

(1.29) 


where a > 1, y is n vector, and f and b are nonlinear mappings 
on suitable domains* Problems on semi- infinite or infinite 
intervals on the other hand usually leads to the case <a 1 
(singularity of the second kind). The author's have introduced 
a canonical form for the system (l*28) which provides a basis 
for the analysis of such boundary value problems* They also 
established a Fredholm theory for linear problems in this 
canonical form and derived existence and regularity results 
for the solutions of nonlinear problems* 


The numerical solution of (l*28) with a = 1, assumed to 
be in canonical form, has been studied by Hoog and Weiss [ '&] * 
Here they extend the results to the general case a 1* The 
centered Euler scheme and the trapezoidal scheme is then used 
to solve the regular boundary value problem* By using Fredholm 
theory, the authors obtain the modified appropriate boundary 
condition for the difference method* They also show for 
important problems the error has an 0(h ) expansion, which 
renders possible the use of Richardson extrapolation or related 
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technique* The appl icati 1 ity of the method has been illustrated 
by Blasius equation and a problem arising in the study of 
electro-magnetism* 

This is the only paper dealing with essential 
singularities that the author has been able to find in the 
. 1 iterature • 

1 *6 Summary of the thesis : 

In the present thesis, some new and efficient numerical 
methods have been developed for solving singular boundary value 
problems over finite and infinite intervals* The thesis comprises 
of eight chapters# 

In the second chapter, a particular class of singular 
boundary value problems have been considered for ordinary 
differential equations on a finite interval, which has a 
singularity of the first kind# These types of problems have 
been discussed earlier by Jamet [50] and arise when reducing 
partial differential equations by physical symmetry* The 
singularity is subtracted out by using the series expansion 
and a new boundary condition is derived in the vicinity of 
the singularity# A discrete invariant imbedding method has been 
developed to solve the resulting regular boundary value problem 
efficiently# The stability of the method is discussed# The 
proposed method is conceptually simple and easy to Implement# 

Some test examples have been solved and the numerical results 
presented# A comparison of numerical results has also been made 
with other methods# 
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Chapter 3 deals with the continuous version of 
Invariant Imbedding method to solve the singular boundary value 
problems considered in Chapter 2* Again the expansion procedure 
has been employed to derive the new boundary condition near to 
singularity* Next a continuous invariant imbedding has been 
described to reduce the regular boundary value problem in to a 
set of initial value problems* Runge Khtta Fehlberg variable 
step size technique has been employed to integrate the stable 
initial value problems efficiently* Some model problems have 
been solved to illustrate the efficiency of the method* The 
continuous invariant imbedding combined with variable step size 
Initial value software have been found to be more accurate 
than the discrete version of invariant imbedding hut at a 
slightly more cost of computational time* 

The treatment given to the singularities to the problems 
considered in earlier chapters is based on using extended power 
series and thereby derive a new boundary condition near to 
singularity* However/ in cases where there is slow convergence 
of the series expansion, an economized expansion may be used* 

In chapter 4, we discuss Chebyshev economization near the 
singular point to overcome slow convergence of the series* 

Using the shifted Chebyshev polynomials and the Lanczos Tau 
method, a new boundary condition near the singularity is 
obtained* The resulting regular boundary value problem is then 
treated with continuous Invariant Imbedding described in Chapter 3* 
Numerical results of some test examples are also included* 
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The attempts In earlier chapters for the removal of 
singularity are based, on using expansion procedures in the 
neighbourhood of the singularity* However for some problems, 
it may be difficult or even not possible bo obtain the series 
expansion near the singularity* In such cases, one may have to 
use a direct approach to solve singular boundary value problems* 
A modified fourth order finite difference method to solve a 
certain class of linear singular BVP / s is presented in Chapter 5* 
The original differential equation is modified at the singular 
point* The main feature of the modified difference scheme is 
that it leads to the tridiagonal system of equations, which has 
been solved by discrete Invariant Imbedding method* Some model 
problems have been solved and the numerical results included* 

The numerical solution of nonlinear singular boundary 
value problems has been discussed in Chapter 6# A quas il inear i- 
zaticn technique has been used to reduce the nonlinear problem 
to a sequence of linear problems* A fourth order finite 
difference method described in Chapter 5 has been employed to 
solve these problems effectively* Some model problems have been 
solved to demonstrate the efficiency of the method* 

In Chapter 7, we have considered the infinite interval 
boundary value problems in ordinary differential equations* In 
Case I, we derive the asymptotic boundary condition to represent 
the infinity condition at the far end* This boundary condition 
has been derived such that it expresses the asymptotic behaviour 
of the solution as well, and converge to the actual solution of 
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the 'infinite' problem as the length of finite interval tends 
to infinity* The resulting boundary value problem is treated 
by imbedding methods* The stability of the problem is analysed* 
In Case II, we have briefly included the transformation technique 
as an alternative method to solve these problems* For the 
transformed problems on a finite interval with a regular 
singularity at one end point of the interval, the method 
discussed in earlier chapters can be applied directly* 

Finally in Chapter 8, the application of cubic spline 
method to infinite interval problems is discussed* By reducing 
the infinite interval to a finite interval which, is large and 
imposing the asymptotic boundary conditions at the far end, 
the two point boundary value problem on a finite interval is 
treated by using cubic spline approximation * The tridiagonal 
system resulting from the spline approximation is efficiently 
solved by the method of imbedding* The stability of the method 
is also discussed and the theory is illustrated by test examples* 

In a nut-shell, the numerical methods presented in 
this thesis for solving singular boundary value problems over 
finite and infinite intervals have been shown to be efficient 
over other methods* Although problems with singularity at one 
end point have been discussed, yet one may apply these methods to 
problems with interior regular singularities or two regular 
singular end points* The invariant imbedding methods is 
conceptually simple, easy to implement, and has the properties 
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of solving boundary value problems subject to different 
boundary conditions without extra computational effort# 

Model problems have been solved and numerical results 
presented in respective chapters « The numerical results 
presented, in this thesis have been computed on DEC-1090 system 
in single precision arithmetic# 
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DISCRETE INVARIANT IMBEDDING FOR SINGULAR 
BOUNDARY VALUE PROBLEM 

2«1 Introductions As mentioned in Chapter 1, the numerical methods 
for approximating the solution of singular two point boundary value 
problems for ordinary differential equations are important in 
view of the frequent occurrence of these problems in Engineering 
and Science. Specific complete algorithms are classified 
according to, how the transformed problem is modeled discretely, 
and how the discrete model is solved efficiently. Due to the 
effect of singularity either in the differential equation or in 
the boundary condition, one may need extra care to solve these 
problems by efficient computational techniques. 

We consider our analysis for the differential equation of 
the form 

Ly = y* * Cx)+p(x)y / (x)+q(x)y(x) = r(x) (2~.l) 

on the interval x Q < x < b, where p(x) and q(x) fail to be 
analytic at x Q , We demand y(x) to be twice differentiable on 
[x Q ,b] and its second derivative to be continuous on (x Q ,b]. 

We also assume that there exists a number cr , 0< cr < 1, such 
that p(x) < — for x small and q(x) < 0, then the two point 
boundary value problem 

Ly(x) = rCx), x Q < x ^ b 

y (x o ) = ° (2.2) 

y(b) = P 

y e C 2 (x 0 ,b] n C[x,b] 
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has a unique solution. For the one point boundary value 
problem 

Ly(x) = r(x), x o < x < k 

yCb) = P (2.3) 

y e C 2 (x Q ,b] O [x Q ,b] 

Id 

with cr > l such that p(x) > — - C , where C is a nonnegative 

**** X 

constant and q(x) < 0, the problem (2.3) has a unique solution. 
For one point boundary value problem we do not impose a value 
on y(x) at x = x Q r the only requir orient is that at the 
singular point y(x) be bounded-; 

A particular class of this singular boundary value problem 
in ordinary differential equation arises as a result of the __ 
application of separation of variables to the equation of the 
generalized axially symmetric potential in a rectangle and 
has been studied by Parter [ 72 ], using finite difference 
scheme. Jamet [ 50 ] studied the problems (2. 2) -(2.3) with 
the application of a standard three-point finite difference 
schane with a uniform mesh size h and has shown certain 
convergence results. Those results appear as particular cases 
of - general results for partial differential equations of 
elliptic-parabolic type which are given by Parter [ 72 ]. 

Gustaffsson [ 47 ] has given a numerical method for 
solving singular boundary value problems by representing 
the solution as series expansions in the vicinity of the , 
singularity and by using difference methods for the rsnain,ing 
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interval. For linear systems of first order, Brabston and 
Keller [ 16 ] have used the sane procedure in the neighbourhood 
of the singularity and used the Box scheme [ 57, 5B ] for the 
solution of regular boundary value problems. 

In this chapter we have presented a computational technique 
for the numerical solution of singular two point boundary 
value problems of the form (2*2). By employing the expansion 
technique on (x Q ,<5], where 6 > x Q , near the singularity, we 
derive a new boundary condition at x = 6* A discrete invariant 
imbedding is then described to solve the problem over the 
reduced interval. The stability analysis of the method is 
discussed. Numerical results for the model problems solved, 
have been presented and their comparisons made with other methods. 

2.2 Removal of Singularity I 

We shall first briefly describe how to replace the 
singular problem by a regular problem on an interval [6,b], 

6 > x Q . We consider a linear singular two point boundary 
value problem 

Ly(x) = y / ' (x)+p(x)y' (x)+q(x)y(x) = 'r(x) ^x Q <x<Cb (2.4) 

subject to boundary conditions 

y(x Q ) = c (2.5) 

y(b) = d 

where x Q is a regular singularity. Under the assumptions 
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pointed out in Sec. 2.1, for the Eqns (2.4)-(2.6) the solution 
exists and it is unique - . 

For the regular singular point x = x Q , we make use of 
the modified series expansion in a small interval near x = x Q , 
so that the Egn. (2 t 4) has a solution of the form (see 
Coddington and Levinson f 2^ ]) 

_ oo 

y(x) = (x-x Q ) p 2 a (x-x Q ) n , a Q f 0 (2.7) 

n=o 1 ' J 

Differentiating (2.7) twice and using in (2.4), and comparing 

the coefficients of powers of (x-x Q ) on both sides, the valuer 

of p as the roots of the indicial equation and the recurrence 

relation for the coefficients a n 's are obtained. Depending 

upon the nature of the roots of the indicial equation, the 

general solution of Eqn. (2.4) can be written as 

n 

y(x) = 2^ a ± S ± (x) + s n+1 (x), n<2 (2.8) 

for x e (x Q , 6], where S^(x) and S^x) are tw ° linearly 
independent solutions of the homogeneous equation corresponding 
to (2.4) and $ n+ ^(x) is the particular solution of (2.4). 

The basic theortical results about these expansions about a 
singular point have been reviewed by Keller [ 5!T ]. The 
series expansion may be valid for the entire interval (x 0 ,b]y 
however due to slow convergence of this series expansion, we 
restrict the expansion in the interval (x Q , §]* where 
x Q < § < x < b. We derive the regular boundary value problem 
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with a new boundary condition as follows. From Eqn. (2*8)/ 
at x = §, 


a^iCg) + a 2 S 2 ( ' 6 ‘ > + S n+1 (6 ^ = Y<6) (2.9) 

and differentiating Eqn. (2.8), at x = 6 * 

ajS^U) + a 2 S^(6) + S^ +1 (6) =^(6) (2.10), 

where the primes denote differentiation with respect to x. 

We solve Eqns. (2. 9) -(2.10) for and oc as 

lY(8)-S n+1 (8Us' 2 ($)-[y' ( / 6)-s' +1 (6)]s 2 (5) 

OC — X - , -- M - « 

Sj (5 ) s * 

_ [y # (6)-S^ 1 (6)]S 1 (6)-[y(5)-S ?l+1 (6 ) ] S|(6 ) 

2 ~ S 1 (6)S'(6) - S 2 (6)S£(5) 


2.11 


( 2 . 12 ) 


Also we have from Eqns. (2.5) and (2.8) at x = x Q , 

a l S l (x o ,+a 2 S 2 (x o ) ^n+l (x o ) = ylx o J 
“l s l Cx o )+<I 2 S 2 (x D ) “ y ( *b' ) - s n+l (3, b ) 


(2.13) 


Using Eqns. (2.11), (2.12) and (2.13), we have 


g(6)s' C6)-g' (6)S 2 (6) g' (6)S 1 (6)--g(6)S£ (6> 

. .... - , v • — S x (x 0 )+ ^ , x ' 


h(6) 


h(6) 


S-(x) 

2 O 


— o-S , (x ) (2*14) 

n+1 o 

where g(x) = y(x)-S n+ ^ (x) and h(x) = S^(x)S' (x)-S 2 (x)S' (x) . 
Eqn. (2.14) can be conveniently written as 


[S' (6)Si(x 0 )-S£(6)S 2 (x 0 )]g(6) + .[S^sJs^a^.) -S^S^CxJ]^ (5) 

= hC5)fc-S n+1 (x 0 )] (2.15) 
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or 

ay(6)+Py*‘ (6) = r (2.16) 

where 

a = [S / 2 (6)S 1 (x 0 )-S£(6)S 2 (x 0 )] (2.17) 

P = ^S^(5)S2(x 0 )-S 2 (6)S^(x 0 )] (2.18) 

and 

7 = hC6)[o-S n+1 (x 0 )]+aS n+1 ( 5 )+ 0S' +1 ( 6 ) (2.19) 

Thus a regular boundary value problem over [_6,b] is given by 
Eqn. (2.4) with boundary conditions (2.16) and (2*63. 

2.3 Discrete Invariant Imbedding l 

The discrete invariant imbedding technique is motivated by the 
fact that it is very simple and easy to implement from the 
computational point of view. Several apparently different 
expositions on invariant imbedding have appeared in the 
literature, (See, for example, Casti and Kalaba [ 19 ], 

Meyer [ 66 ], Angel and Bellman [7 ]). 

We describe the method of discrete invariant imbedding 
for the reduced regular boundary value problem given by 
Eqng,(2.4), (2*16) and (2.6). A unique solution of the problem 
given by ( 2 . 4 ), (2.16) and (2.6) exists if p(x), q(x),r(x) 
are continuous on [6,b] and q(x) is negative there. Since 
these functions are continuous on a closed bounded interval, 
there must exist positive constants p*, cf* and q^ such that 

lp(x)I p*, 0 <(^<.q(x) < q* , 5 <x <b. 
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We now divide the interval [§,b] into n equal parts as 
d ~ x Q ^ ^1 ^ x 2 * * * ^ x n = b where x^**x^ ™ h, (r = l,2,.*,n). 

Employing central difference approximations for the first and 
second order derivatives, the Eqn. (2.4) is discretized as 


^i+l-^i^i-l 


+ p. 


Y i + l~ Y i-l 


2h 


+ q i Y i = r i 


i =0,1,2,.. ,n-l 

( 2 . 20 ) 

where p. /q. and r. denote p(x. ), q(x. ) and r(x. ) respectively. 

The difference Eqn. (2.20) can be written in a tridiagonal 
form as. 


~ A ±yi+i + B ± y i - c i y i-i = D i 
where 


A 


and 


± = 1 + hp ± /2 


B ± = 2 - h"q ± 


C i = 1 - hp ± /2 


D i = h r i 


x ss 0, 1, 2, • « , n— 1 (2.2l) 

( 2 . 22 ) 
(2-. 23) 
(2„24) 


(2-. 25) 


Consider now the coupled difference equations by writing 
Eqn. (2.21) as 


v i + l * Y i 


(2.26) 


A+i = M ± Y i - N i v i -°i 

where = B.^ ,* N ± = arfi 0. = D../A., 


(2.27) 
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The boundary condition are discretized as 


(Yi-y_i) 

ay o + P - 2h- - - r 


(2.27) 


Y n = d 


We seek a solution to Eqns. (2. 25 ) - (2 . 26) in the form 


v i - W i y i +T i - 


i = 0, 1, 2 , . .,n-l 


(2.28) 


(2.29) 


where W^ and T^ represents w(x^) and T(x^) respectively. 

By rewriting Eqn. (2.29) and using (2.25)*-(2.26) to relate 
W ±+1 and T^ + ^ with and T^ we have, 

n = w i+i <M i y i - N i v i-°i ) + T i + i 

= « i+ i (M i y i - N i (w i y 1 +T i > -°i > « i+1 

= (W ' _ (M.y.-W. ,-N.W.y. ) + (T. , , -W. . 0 . -W. ,.N . T. ) 

r+1 ii l+l i ri i+l i+l i i+I i i 

= (w. , -M.-W. , ,N.W. )y. + (T. L - -W. 0 . -W. ,.N.T.) 

l+l l l+l l i l+l i+l i i+l i i 

we match the coefficients in y^ to obtain the relations 


i+l w. ) 

1 11 


(2.30) 


T. _ = (W.,,0. + W. ,.TJ.) 

i+l l+l i i+l l i 

The Eqn. (2,30)-(2.3l) give the recursion relations 


(2.31) 


for W ±+1 and T ±+1 for i = 0,..,n-l. To obtain W i+1 and ^ ^ 
we need to know the value of W^ and T^ at i = 0. To do this 


we have from Eqn. (2.27) as, 

2hr , 2 ha 

y i - v + y _i - y o 


(2.32) 
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But from Eqn. (2,21) at x = x Q ( — &) we have 


“ A o Y l + B o Y o “ C o Y -l = D o 
Now using (2,32) and (2.33) we get. 


(2-. 33) 


- A o[ 

2 hr , 

nr + 

,, 2 ha 

Y -1 “ P* 

Y ol 

+ 

By 

0^0 


2A ha 





Y o[ 

o 

'P " 

+ B oH A 0 +c ofc 

-1 

= [ D o + 



2A ha 



Y -1 

t A o +C o] 

= [B q + ■- 

O 

P ' 

bo 

.-cv : 


2A hr 


2A hr 

o 


2A ha 

L V ■ -3 


y o - 


2A hr 

[ D o+- pP- 


[A +C ] [A +C ] 

u o o J L o o J 

Comparing the Eqns. (2.29) and (2.34), the initial values 


(2.34) 


of W and T are obtained as, 

2A ha 

tv-'M 

w ° = 


(2.35) 


T = 
o 


2A hr 

-LV 


( 2t36 ) 


From the Eqns .(2. 35 )- (2, 36 ) and (2„30),(2.3l) we obtain the 
values of VA and in the forward process from 

i = 0, 1, 2, . . ,n-l. To obtain the solution y^ , we have 
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y. ,, = M.y.-N. v.-O. 

J i +1 i J i 11 i 

= M.y.-N. (W.y.+T. )-0. 
x X X 1 x 


= (M.-N.W. )y.-N.T.-0. 
1 11 -*1 XX X 


y± = 


(y. , -+0.+T.N. ) 
fi+X i i_j. 

(M.-N.W. ) 
l 11 


(2.37) 


This equation can be used to find y^ where i ranges from 
n-1 to 0, in the backward process using the stored values 
of and T^ and the boundary value y n = d. 

2.4 Computational Algorithm l 

In order to compute the solution y^' s (i = 0 / ..,n-l), the 
computation is performed sequentially as follows t 

Step 1 l Find the values of a, P and 7 from Eqn. (2.19). 

Step 2 l Calculate the values of W Q and T q from Eqns. (2.35)- 

(2.36) and find the values of W^ + ^ and 
(i = l/2/../n-l) in the forward process using the 
Eqns. (2.30)-(2.3l) 

Step 3 ’• Compute the value of y/ s (i = n-l / .. / 0) by • 
using Eqn'. (2.37) and the known W. 7 s and T. / s j 
(i = 0, ,,.,n) from Step 2 and using Eqn. (2.6),. 


2.5 Stability l 

% 

By stability we mean that the error made in one stage 
of the computations is not propagated into larger error 
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in succeeding calculations. Now we examine the stability 
for the recurrence relations given by Eqns. (2.3°) and (2»3l). 


Let us assume E^ be a small error that has been introduced 


in the calculation of VL , we then have 


W jL = W i +E ± 


(2.38) 


Thus instead of using Eqn. (2.30), we actually solve 


W 


i+1 


[Hi-N-Wj 


(2.39) 


From Eqns. (2 _ . 38) and (2.39) we find the error at the 
next stage is given by 

E i+i = - va 1 ' 1 

= [M i -N i (w i +E i )]~ 1 -(M ± -N i W i )*' 1 
= (M.-N.W.-N.E. ) _1 -(m.-n.w. ) _1 

1 1 1 1 X XXX 

= (M.-N.W.-N.E. )“’ 1 N.E. (M.-N.W.) -1 

111X1 1X1X1 

~ (M.-N.w.) -1 N.E. (M.-N.W.)” 1 

“■ill lllll 

.2 


= W~ C./A. E. (2.40) 

l+l x / l 1 

under the assumption that the error is initially small. 

Assume that N^> 0 for O < i < n-1 and since q(x) < 0 it 
can be easily verified that > A^+Ch , for 0<i _< n-1. From 
Eqn. (2.35) we have 

9 A hY 

W o - < B 0 + ■ g- ] 
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2 

and 1 W Q 1 < 1 if M > 2, where M = h 2 ^- - ^ r p o (2.41) 

under this condition and making use of the assumption on 
N i / it follows very easily from (2.30) that 

iWj <1, for i = 1,2, ...,n. (2.42) 

From Eqn, (2.40) it then follows that 

IE^ +1 1 < lE^l , provided !cJ < 1 A J 

and thus making the recurrence relation stable. 


Similarly for the other recurrence relation we have 


T i - T ± + e ± 


where is a small error, from Eqn. (2.3l), 

*i+i = w i + i°i%iVi 


(2.43) 


(2-. 44) 


We find the error in the next stage is given by 


£ i+l " T i+1 ~ T i+l 


= W. . e .N. 
i+l i x 


and so 


[ a I < le^i # IW i l < 1, for i = 1,2, ».,n and the 
assumption |C^l<lA^i making the relation (2.31) stable. 

2.6 Numerical Results and Discussion ! 

In this section, the numerical results for the model 
problems are given. All computations were carried out in 
single precision arithmetic on DEC- 1090 computer systen. 
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Problem 2.1 l We consider the linear two point boundary 
value problen 


y" (x) + ^ y' (x)--cyCx) = 0 , 
with boundary conditions 


O < cr< 1 
r >0 


y(o) = 1 
y(l) = 0 


This is the model problem solved earlier by Jamet [ 50 ]. 


Problem 2.2 \ Next we make the numerical experiments on 
the equation 

y // (x) + — y' (x) = -x 1 " 0 cosx-(2-o)x"~ 6 sinx, 0 < o < l 
with boundary condition 


y (0) = 0 

y(l) = cosl. 

This is the model problem considered by Gustaffsson [ 47 
which has the exact solution y(x) = x~~ a cosx. 

Problem 2.3 * Consider the linear two point boundary value 
probl em 

(xV)' = Px a+P " 2 ((a+P>l) + P x P )y , 0 < x <1, 

with boundary conditions 

y(0) =1 
y(l) = e. 
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This model problem has been taken from Chawla [ 21 ~] f which 

D 

has the exact solution y(x) = exp(x P ). 

Problem 2.4 * We consider the problem 

(x 1/2 y' )' - x 1/2 y = - |(3+x" 1/2 )e" X / 0 < x< 1 

subject to boundary conditions 

y(0) = 1 
y(l) = 2e~ 1 . 

This is the model problem solved by Doedel and Reddien [ 37 ] / 
which has the exact solution y(x) = (x^ 2 +l) e'” x . 

Problem 2.5 As a last model problem, we have 

Y" - y' = -1 , 0 < x <1 

X '■*** 

subject to boundary conditions 

y(0) = 5 
,y(l) = 5 

9 * 5/0 

whose solution is y(x) = 5-x +x ' . This model problem 

has earlier been discussed by Ewa Weinmuller [92,93]« 

- The numerical results 'for problem 2.1 at "different mesh 
points for several mash sizes and two different values of 
6 are presented in Tables 2.1 -2.2 . It can be found 
that these results are much accurate than the results of 
Garnet [ 50 ] . To make it more apparent the comparison of 
our solutions for different mesh sizes is made with two 
other methods and presented in Table 2.5. It can be 
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observed that the numerical results obtained by the present 
method are quite comparable for any mesh size with the finite 
difference method of Jamet [ 50 ] and the collocation method 
of Reddien £ 73 ]• The solution at x =0.5 for cr =0.5 correct 
up to five significant places is known to be 0.25203 for this 
problen. Tables 2,3 - 2.4 present the solutions for the same 
example for a different value of o and two choices of values 
of 6, For this value - of o t the comparison of the solution 
is made with the Jamet' s solution [ 50 J and presented in 
table 2.6. 

In Table (2,7), the comparison of maximum errors in the. 
solutions for problen 2.2 for different values of 6 is made, and 
the errors in Gustaffsson' s solutions correspond to the finite 
difference scheme of order two. The superiority of the solutions 
obtained by our method is again evident from the results. 

In Table 2.8 the numerical results for the problem 2,3 
are presented for different values of h» The results are 
given for P = 4 and a = 0.5. To test the effect of singularity 
.the values of 6 are taken near to singularity. For problems 2.4 
and 2.5, the results are shown in Tables 2*9 and 2.10. It, can 
be observed from these tables that the computed solutions , 
compare "well with the exact solutions. 
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Table 2.1 

Results for Problem 2.1 




(cr = 0.5, 

H 

II 

H 

•* 

O 

o> 

= 0.5) 


■ 


h 

- ■ 


-W-k. 

X 

1/10 

1/20 

/ 

1/40 

1/80 

1/160 

0.5 

0.2509616 

0.2517711 

0.2519739 

0.2520425 

0. 2520445 

0.7 

0.1363321 

0.13 67357 

0.13 68365 

0.1368616 

0.1368656 

0.8 

0.0873957 

0.0876477 

0.0877106 

0.0877263 

0.0877285 

0.9 

0.0422726 

0.04 23 922 

0.0424220 

0.0424294 

0.0424304 


Table 2.2 

Results for Problem 2.1 
C a = 0.5, t = 1.0,6 = 0.2) 

h 


X 

, 1/10, 

. ..1/20 ; . 

1/40 

1/80 _ 

1/160 , 

0.2 

0.5019590 

0-. 5065314 

0.5076795 

0.5079669 

0.5080301 

0.5 

0.2495478 

0.251^216 

0.2518862 

0.2520021 

0.2*20220 

0.6 

0.1892235 

0. 1906261 

0.1909706 

0.1910565 

0.1910704 

0.8 

0.0869033 

0.0875260 

0.0876801 

0.0877185 

0.0877241 

0.9 

0.0420345 

0.0423333 

0.0424Q72 

0.0424256 

0.0424282 
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T eble 2. 3 

Result s fo r Problem 2.3. 




(Cf=0.2 , 

t= 1.0/ 6 = : 

3.5) 


X 



h 




1/10 

1/20 

1/40 

1/80 

1/160 

0. 5 

0.3727644 

0.3732550 

0.3733774 

0.3734079 

0.3734118 

0.7 

0.2117490 

0. 2120046 

0. 2120681 

0.2120838 

0.2120841 

CO 

• 

o 

0. 1383073 

0.1384698 

0.1385101 

0.1385201 

0.1385198 

0.9 

0.0680450 

0.0681234 

0.0681428 

0.0681476 

0.0681473 


Table.. 2.4 

Results for problem 2.1 




( a = 0.2, 

T=1.0/ 6 =0 

.2) 


X 



h 



1/10 

1/20 

1/40 

1/80 

1/160 

0. 2 

0.6756758 

0.6770162 

0.6773548 

0.6774390 

0.6774540 

0.4 

0.4632033 

0.4639039 

0.4640781 

0.4641208 

0.4641227 

0. 5 

0.3727149 

0.3732428 

0.3733736 

0.3734055 

0.3734049 

0.6 

0. 2894291 

0.2898205 

0. 2899173 

0.2899407 

0.2899389 

o 

• 

00 

0.1382890 

0.1384653 

0.1385087 

0.1385192 

0.1385174 

0.9 

0.0680360 

0.0681212 

0.0681422 

0.0681472 

0.0681462 
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Tabl e 2._5 




Comparison of 
(<7= 0.5, 6 = c 

Results for the problem 2.1 

1 . 5, T= 1) 

N = 

1 

h 

J arnet' s 
Method [50 ] 

Reddien' s 
Method [73 ] 

Discrete Invariant 
imbedding 

8 


0. 2903821 

0.25305 

0. 2503571 

16 


0. 2782581 

0.25223 

0.2516190 

32 


0.2700966 

- 

0.2519360 

64 


0.2645623 

- 

0.2520172 

128 


0. 2607722 

- 

0.2520350 

256 


0.2581542 

- 

0.2520371 



Table J?« 6 




Comparison of 

Results for the problem 2.1 


( <7 = 0.2, 6 = 0.5,t = 1.0) 


N = Discrete Jamet's 


h 

Invariant imbedding 

Method [50] 


8 

0.3723980 

0.3830 297 


16 

' 0.3731628 

v 0.3789884 


32 

0.3733548 

0.3766250 


64 

0.3734027 

0.3752605 


128 

0.3734142 

0.3744761 


256 

0.3734141 

0.3740257 
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CHAPTER hi 


CONTINUOUS INVARIANT IMBEDDING FOR 
SINGULAR BOUNDARY VALUE PROBLEMS 

3*1 Introduction: The use of invariant imbedding techniques* in 
general* reduces the linear two point boundary value problems 
into a set of initial value problems* In the earlier chapter* 
a convenient discretized version of invariant imbedding has 
been developed for solving singular boundary value problem after 
its reduction to a regular problem by deriving a. new boundary 
condition in the neighbourhood of the singular point* Though 
a reasonably good accuracy is obtained for a moderate step size 
h, yet one needs to take* as expected, a very fine mesh to 
achieve high precision in the solution especially when the point 
x = 6 is very near to the singularity* This is due to the 
fact that the difference scheme of 0(h ) would not, in general, 
give solutions correct upto second order near the singular 
point* The choice of a small step size is undesirable due to 
the fact that it leads to computational complexity* 

To overcome these difficulties, we present another 
computational method on the basis of expansion technique and 
continuous version of invariant imbedding for solving singular 
boundary value problem considered in Chapter 2* The application 
of continuous version of invariant imbedding has been discussed 
in the literature by many authors* (see, for example. Bellman and 
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Kalaba [l2], Chandrasekar [ 20 ], Bellman and Wing [13], Nelson 
and Scott [69], Kalaba and Kagiwada [53] and Denman and Mingle [34] 

Sufficient conditions for existence of two versions of 
invariant imbedding for linear second order equations with a 
regular singular point have been considered by Scott [85] • The 
work of Scott constitutes a significant extension of earlier 
results of Banks and Kurowski [lo]» Elder [38] showed how to 
apply the subscheme of invariant imbedding known as integration— 
to— blowup to compute the smallest eigenlength of linear first 
order system having a singularity of first kind* Nelson, Sagong 
and Elder, [7l] adapted Elder's approach to the solution of certain 
singular two point boundary value problems by means of Scott's 
version of invariant imbedding* Specifically, this approach 
applies to problems based on a linear homogeneous first order 
differential system with a singularity of first kind and with 
boundary conditions consisting of existence ( finite) at the 
singularity and specified values at some second point* 

In this chapter we have presented a numerical method 
for solving singular boundary value problem with a regular 
singular point at one end of the interval* The method is based 
upon the series expansion technique and continuous invariant 
imbedding* A new boundary condition at x = 6 is derived by 
using series solution in the vicinity of singular point to 
substract out the singularity* The regular boundary value problem 
is then solved by employing a continuous invariant imbedding 
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method* Some numerical experiments have been included to 
demonstrate the applicability of the method* 

3*2 Reduction to Regular Problem : 

We shall consider a second order linear singular 
boundary value problem of first bind given by 

y"(x} + F 1 (x)y / (x)+P 2 (x)y(x) = F 3 (x), < x < b (3*l) 

subject to boundary conditions 

y(x^) * a (3*2) 

y(bj = 3 (3*3) 

where a and 3 are given constants* 

Since x = Xp is a regular singular point, we make use 
of the modified series expansion in a small interval (x^,6j so 
that the Eqn* (3*1) has a solution of the form 

CO 

y(x) = (x-x^) P S a n (x-x Q ) n / a Q / 0 (3*4) 

n=0 

Differentiating (3*4) and substituting in (3*1) and comparing 
the coefficients of powers of (x~x Q ) on both sides, the values 
of p as the roots of the indie ial equation and the recurrence 
relations for the coefficients a^/ s are obtained* Depending 
upon the nature of the roots of the indicial equation, the 
general solution of (3*l) can be written as 



56 


m 

y(x) = R ± (x) + R^tx), m < 2 ( 3 * 5 ) 

for x e (3^/dJ , where R^(x) and R 2 (x) are two linearly 
independent solutions of homogeneous equation for ( 3 *l) and 
R m+1^ is the Particular solution to (3.1)* By taking the 
series expansion on (x^/ 6 j# the new boundary condition at 5 is 
derived* We have from Eqn* (3*5), 

a l + a 2 R 2 ( ' <5 ^ + R m+1^ <5 ^ = Y ^ 6 ) (3*6) 

a l + a 2 R 2 (5) + R m+ 1 1 2 * * * (6) = Y<(5) (3*7) 

where primes denote the differentiation* Equations ( 3 * 6 ) — (3*7) 
are solved for and a 2 as 

1 ) -'R2<« 

and 

Cy'Ce )-R 4 fX Ca )] ^(6 )-[y (6 ) ~Vi (6 ^ R i ( ?l 

2 R^d ) R^OS) - R 2 ( 6 ) R£( 6 ) 

Also from Eqns* (3*3) and (3*5), we have 

a i W + a 2 R 2 (x b ) + = y(Xo> (3#10) 

Using the equations (3*8), (3*9) and (3*10) and using 

y(x c ) = a, and after simple algebra, we have 


( 3*8) 


(3.9) 



where 


A y(6 ) + B y 7 ( 6 ) = c 


(3-11) 


A “ t R 2^ 6 ^ - Ej(6) ^(x^,)] 

B = [^(6) Rj^) - Rj(6) HjC^)] 

C = h(6) [»-y % )] + A Vl (5)+B R m + 1 (6) (3 ’ 1S) 

with 

h(6) = [Rj_(6) RJ(6) - 1^(6) R£(6 )] 

The Eqns* ( 3 •! ) (3*ll) and (3*3) form a reduced regular boundary 

value problem over the interval [6,b]* 

3*3 Continuous Invariant Imbedding : 

In this section/ we discuss the method of continuous 
invariant imbedding for solving the reduced problem over [6/b] * 
We follow on the lines of Scott / s version of Imbedding [85,87] 
to reduce the boundary value problem into a sequence of initial 
value problems* 

To be specific, we write Eqn* (3*1) as a systems of 
two first order equation in the form y(x) = u(x), (y(x)) = 

(u(x)) = v(x) and — = ~T^(x)v(x5-F 2 (x)u(x)+P 3 (x)* Hence 
we may write Eqn* (3»l) in the form 

duJLsA = v(x) (3*13) 

dx 

- = P 1 (x)v(x)+F2(x)u(x)~P 3 (x) ( 3*14) 

where the minus sign on the left hand side of the Eqn. (3*14) is 
purely historical# The corresponding hoandary conditions are 
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A u(6) + B v( 6 ) a C (3.15) 

u(b) a 13 (3.16) 

We now discuss the invariant imbedding for a more general’ first 
order system of the form 

u 7 (x) = a(x)u(x)+b(x)v(x)+S + (x) (3*17) 

~v 7 (x) = c(x)u(x) + d(x)v(x) + S~'( x) (3*18) 

where primes denote differentiation/ with all functions 
a(x) , b(x)/ c( x) / d(x)/ S + (x) and S“(x) to be continuous and 
subject to boundary conditions given by (3*15) - (3*16)* From 
the method of continuous invariant imbedding/ we have the 
relations (see Scott [35])* 

u(x) a S 1 (x) v(x) + S 2 (x) u(fi) + S 3 (x) (3*19) 

x e [fi/b] 

v(6 ) = Q^(x) v(x) + Q 2 (x) u(6 ) tQ^Cx) (3*20) 

where the Eqn* (3*19) is the generalized Ricatti transformation 
and Eqn. (3.20) is the recovery transformation. The first term 
on the right hand side of Eqn.(3*l) represents the contributions 
from the homogeneous equation/ the second term represents the 
contributions from the boundary conditions at x = 6 and the 
last term represents contributions from the nonhomogeneous term 
in the equations* The transformations can be motivated from both 
a physical and mathematical point of view* 



We shall now derive differential equations for the 
functions S^x), S 2 (x), SgCx), Q^x), Q 2 (x) and q 3 ( x >* 
Differentiating (3-19) we have 

u'(x) = S'(x)v(x)+S 1 (x)v / (x)+S'(x)u(6 )+S'(x) C 3*21 ) 

using Eqns. (3-17) and (3-18) in (3*2l) we obtain, 

a( x)u(x)+b(x)v(x)+S*(x) = S^(x)v(x)+S^(x){~c(x)u(x) 

~d(x)v(x)-S (x)} + Sp( x)u(s )+S^( x) 
or 

a(x){S^(x)v(x)+S 2 (x)u(6)+S 3 (x)} + b( x)v( x)+S + ( x) 

= s£( x)v( x)+S^( x) [ -c(x){ S^(x)v(x)+S 2 (x)u(6 )+S 3 ( x)} 

— d( x)v( x)—S (x)] + S 2 (x)u( 6) + S^(x)« (3-22) 

or 

v( x) { S^( x) a( x) +b( x)~S£( x) +c( x)S^( x)+d( x)S^ ( x) } 

+ u(d) { S 2 ( x) a(x)+c( x)S^( x) S 2 ( x)-S 2 ( x) } (3-23) 

+{S 3 (x)a(x)+S + (x)+S 1 (x)S 3 (x)c(x)-fS~(x)s 1 (x)~S^(x)} » P 

Eqn- (3-23) will be satisfied if each term in the braces to set 
equal to zero* That is 

2 

S'(x) = b(x) + [a(x)+d(x)] S^(x)+c(x)S^( x) 


(3*24) 


DU 


S'(x) = [a( x ) + c(x) S 1 ( x )] s 2 (x) (3.25) 

S'(x) = [a(x)+c(x)S 1 (x)] S 3 (x)+s‘~(x)S 1 (x)+S + (x) (3-26) 

suitable initial conditions suggested by (3.15) and (^.16) are 


S^Cfi) = 0 f S 2 (6) =i ; S 3 (6) = 0* (3-27) 

Similarly differentiating (3*20) and using Eqns- (3-17) ~ (3*18), 
we obtain another set of equations/ 

Q^(x) = [d(x) + c(x) S 1 (x)] Q^x) (3.28) 

Q 2 (x) = [c(x) Q^(x) S 2 (x)] Q 2 (x) (3*29) 

Q'(x) = [ c( x) S 3 (x) + S~(x)]Q 1 (x) (3-30) 

with suitable initial conditions 

Q 1 (6) = 1 ; Q 2 (6) =0 ; Q 3 (6) = 0 (3.31) 


One can easily verify that the above process does indeed solve 
the original linear problem. Thus sequence of initial value 
problems (3*24) - (3*31) are solved numerically by efficient 
initial value routines* Once S^(x), S 2 (x), S 3 (x), Q^(x), Q 2 (x) 
and Q 3 (x), x e [6, b] are known, the Eqns. (3*19) - (3*20) can 
be evaluated at x = b to yield 

u(b) = S^(b) v(b) + S 2 (b) u(6) + S + (b) (3*32) 

v(6) =Q 1 (b) v(b) + Q 2 (b) u(«) + S~(b) (3-33) 
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Since u(b) is known from Eqn, (3*16), we solve Eqns. (3-32), 
(3*33) and (3*15) as a system of linear equation for the 
unknowns u(6 ) , v(6 ) and v(b). Then Eqns- (3*19) and (3*20) 
are solved for u(x) and v(x) for all values of x £ [d,b], by 
using the known values of u(fi), v(s) and v(b) respectively. 

3 « 4 Computational Algorithm : 

In order to compute the solution u, the computation is 
performed sequentially as under : 

Step 1 : Integrate the initial value problems given by Eqns* 

(3*24) - (3*31), using efficient initial-value routines, 
from x = 6 tox = bto obtain the S^ and profiles 
and store them. 

Step 2 : Find the values of A, B and C from Eqn. (3*15). 

Step 3 : Evaluate the Eqns. (3*19) — (3*20) at x = b using the 
values s j_(b) and Q^(b), (i = 1,2,3) from Step 1, and 
using u(b) from the boundary conditions* The resulting 
equations would contain the unknowns u(d), v(fi) and 
v(b) * 

Step 4 t Combine the equations resulting from Step 3 with 

equation .( 3 *15 ) and solve these three equations for 
the unknowns u(6), v(s) and v(b). 

Step 5 s Compute v(x) and u(x) for any x £ [fi,b] using the 
values of u(6) and v(6 ) from Step 4, and the stored 
values of S^ and from Step 1 • 
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It may be noted that the Step 1, which incidently involves 
maximum computer time in the entire computational procedure/ 
need not to be repeated if one has to solve different boundary 
value problems given by the same differential Eqn-(3*l) but 
different boundary conditions (3*2) - (3*3) but only Steps (2) 
(6) be repeated for each set of boundary conditions, using the 
stored values of and from Step 1 in each case* 

3 *5 Numerical Results and Discussion t 


In this section we shall illustrate the use of the 
algorithm derived in the previous section by applying it to 
several examples- In our examples all the corresponding initial 
value problems (^#24) - (3*2l) were solved by using a fourth- 
fifth order Runge Kutta Fehlberg scheme designed to estimate the 
local error and control step size for the accuracy requirements 
[40,41] * All computations were carried out in single precision 

on DEC-1090 computer system with a relative and absolute error 

-5 -13 

of 10 and 10 respectively - 

Problem 3*1 : We consider the linear two point boundary value 
problem 


u"(x) + - u'(x) - Tu(x) =0, 0 < a < 1, 

A T > 0 


with boundary conditions 


u(0) = 1 


u(l) = 0 
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This proolem has been studied earlier by Jamet [ 50 ]* 

Problem 3*2 : We solve 

2x(l-hx) y* + ( l+5x) y** + y = 0 
with boundary conditions 

y(Q) = 1 
y(l*5) = l 

This problem has earlier been solved by Cohen and Jones [26] 

and has the exact solution given by y(x) = . 

( 1 + x) 

Problem 3*3 : Next we make numerical experiments on the non- 
homo gen eous equation 

rx A | i , L 0 

u'^x) + — u'(x) = -x x ~ cos x-(2-cr)x sin x* 0<a <i 

with boundary condition 

u( 0 ) =0 
u( l) = cos 1 » 

This example has been taken from Gu staffs son [47] and has the 

. * l-a 

exact solution u(x) = x cos x* 

The computed solution for Problem 3*1 for different 
values of 6 when o = 0*5 have been presented in Table 3*1* At 
the various values of x* the solutions compare very well with the 
analytical solution (series solution)* The problem (3*1) has 
been earlier solved by Jamet [50] and Reddien [73] and their 
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computed solution at x = 0*5 for different mesh sizes have 
been given in Table 3*2 •’ This table also contain solutions 
computed at x = 0 *5 by Invariant imbedding using a Runge—Kutta— 

F ehlberg scheme with step size control to solve our initial 
value problems- Thus the step sizes h given in the Table 3*2 
have no direct relevance to invariant imbedding solution* As 
in evident from this table, the invariant imbedding solution is 
much superior to the solution obtained by Jamet or Reddien 
for a very fine mesh- The behaviour of the solution when 6 is 
very near to the singularity has been shown in Tables (3-3). ~ 
(3*4)* It has been observed here that a comparatively smaller 
step size is required (as expected) to achieve the desired 
accuracy in the solutions* 

The numerical solutions for problem 3*2 obtained by 
our method and the comparison of those solutions with Cohen's 
solutions have been given in Table (3-5)- The solutions 
obtained by our method compare well with exact solution and that 
of Cohen's solutions* 

Tables (3*6) - (3*7) give the numerical results for the 
problem 3*3* The computed solution at different points -with 
respect to several values of 6 are presented in Table 3*6* The 
comparison of maximum error incurred in the solutions obtained 
by our method and that of Gustaffsson has been made in Table (3*6)* 
The superiority of our method in most of the cases is evident 
from the results * 
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TABLE 3.1 

Results for the Problem 3.1 

<5 = 0-2 6=0. 4 6=0. 5 

0.2 0.5030600 

0*4 0.3221206 0.3221289 - 

0 0.25,20316 0.2520420 0.2520425 

0.9 0 .0424309 0.0424321 0.0424322 


TABLE 3» 2 

Comparison of Numerical Results for Problem 3*1 
(o' = 0*5, T = 1, X = 0.5) 


N = 

1 

K 

Jamet' s 
Method [50] 

Reddien' s 
Method [73] 

Invariant 

imbedding 

Method* 

8 


0.29038211 

0.25305 

0.25 204250 

16 


0*278 25809 

0.25223 

— 

3 2 


0.27009658 

- 

— 

128 


0.26077 219 

- 

~ 

512 


0.256 33371 

— 

— 


*The step size h has no direct relevance to Invariant imbedding 
solution and refers only to Jamet^s and Reddien' s solution. 



TABLE _3 • 3 

Results for the Problem 3.1 
(6 = 0.05, a = 0*5, T = i) 


X 

u(x) 

0 .05 

0.75006408 

o 

♦ ■ 

o 

0.64846210 

0.15 

0.57142780 

0.20 

0.50806241 

0.50 

0.25204202 

0 .80 

0.0877 3196 

0.90 

0.04243212 


TABLE _3 .4 

Results for the Problem 3.1 
( 6 = 0.125, a = 0.5, T *r 1) 



u(X> 


0.125 
0.250 
0.375 
0 .500 
0.6 25 
0 .750 


0.60761680 
0 .45 34 39 30 
0-34147399 
0.25204244 
0.17697532 
0.11175222 


0.875 


0.53454321 



67 


TABLE 3.5 



Results 

for 

the 

Problem 3-2 ( 6 = 0.5) 


X 

Imvariant 

Imbedding 

Solution 



Cohen' s 

Solution 

Exact 

Solution 


... , _ 




,T ■ >. - ' if. . ' ' . ■ W 

0.5 

1.243902 



1 -2439 6 4 

1*244017 

0.6 

1.217836 



1.217871 

1*2179 21 

0.7 

1.1909 24 



1.1909 39 

1-190997 

0 .8 

1.64078 2 



1-164078 

1.164136 

0.9 

1.13779 4 



1.137781 

1 *137840 

1*0 

1-112338 



1.112319 

1.112372 

1 .1 

1.087342 



1.0878 24 

1.087368 

1 *2 

1.064364 



1.064348 

1.06438 2 

1 .3 

1.041912 



1 .041900 

1.0419 24 

1 .4 

1.020469 



1 .0 2046 2 

1.020474 




TABLE 3-6 



Results 

for 

the 

Problem 3*3 (a =0.5) 


X 

6 = 0 .1 



6 = 0*2 

6 =0.4 

0.1 

0.31464841 



- 

- 

0.2 

0.438 29967 



0 .438 293 9 2 


0.3 

0.52326032 



0.52325916 

- 

0 .4 

0.58 253142 



0*58252991 

0 .58 25 301 

0.6 

0.639 30355 



0.6 39 30194 

0.6 39 30 2 2 

0 .8 

0.62315431 



0.62315313 

0.62315 34 

0 .9 

0.53971162 



0.58971087 

0.5897109 
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TABLE 3.7 





Comparison of Maximum 

Errors* 

N = 

1 

K 

0 

Gustaffsson 

Invariant Imbedding 




Solution 

Solution 

40 



1 »5xlO~ 5 


80 


0.1 

7 «2xlO*~ 7 

4«4xl0 8 

160 



4*OxlO*~ 8 


40 



7 •9xlO~ 7 


30 


0.2 

4-4xlO~ 8 

1.8xl0~ 8 

160 



2-6xlO~ 9 


40 



3-7xlO~ 3 


80 


0.4 

2.1xl0“ 9 

7 -4x1 o" 9 

160 



1 •3xl0*' 1 ° 



*The step size h refers only to the Gustaffsson case* 



CHAPTER IV 


Singular boundary value problems Via Chebysbev 
Polynomials and Invariant imbedding 

4'«i Introduction; The object of this chapter is to describe 
a method based on Chebyshev polynomials and Invariant imbedding 
to singular boundary value problems given in Chapter 2* We 
have observed that for a particular class of problems the 
series expansions may be suitable due to the nature of the. 
physical situation of the defined problem., but in general it 
alone may not produce an effective approximation near the 
singularity# It is due to the fact that the series may 
converge slowly in the vicinity of the singularity thereby 
necessitating to include a large number of terms of the 
expansion to obtain a reasonably good accuracy* However an 
attempt is made here to reduce the necessary number of terms 
in the series expansion and without increasing the errors by 
the process called Chebyshev economization* This method 
effectively produces a good approximation and is based on the 
fundamental property of the Chebyshev polynomials » 

Recently Cohen and Jones [26] studied a shifted 
Chebyshev polynomial with finite difference correction approach 
for a second order linear ordinary differential equation with 
a regular singular point* They considered these polynomials 
on the whole interval where the polynomials are valid, by 
neglecting the effect of singularity* 



In this chapter we study the economized expansion 
procedure based on Chebysnev polynomials to remove the 
singularity for a singular two point boundary value problem* 

The new boundary condition is derived in the vicinity of the 
singularity* The invariant imbedding method discussed in 
Chapter 3 is used for numerically solving the reduced boundary 
value problem* Model problems have been solved and the 
numerical results are presented. 

4*2 Economized Series E:xpansion : 

For the sake of brevity, we consider a homogeneous 
linear ordinary differential equation given by 

Ly = y 7/ ( x) +F^( x) y* ( x)+F 2 (x)y( x) =0, 0 < x < 1, (4*l) 

with boundary conditions 

y(0) = a (4*2) 

yll) =3 ( 4 * 3) 

where x = 0 is the regular singular point of the differential 
equation ( 4 *1 ) * However if the interval under consideration is 
[a,b ], we transform it to the interval [—1,1] by the change o±. 
variable, 

x = ^ (b-a) t + -^ (b-f-a) (4*4) 

and for the special range 0 < x < 1, we can write 


t = 2x - 1 


(4*5) 



Due to the singularity at x = 0, the solution of Eqn.(4.1) can 
be written as 

OO 

y(x) = x r £ C x n , C fc 0 (4 »6) 

n=0 ° 

The coefficients C^'s and the indicial roots ris can be obtained 
by differentiating (4-6) and substituting in (4 -I), comparing 
the coefficients of powers of x on both sides of the equation* 

One may write the general solution as 

2 

y( x) = s a ± R ± (x) (4*7) 

i=l 

where R^(x) and R^tx) are two linearly independent solutions 
and and are arbitrary constants* In a situation where 
we have slow convergence of the series solution,, we express the 
solution of (4-1) by a Chebyshev Economized expansion in (0,s] , 
where s is chosen near to singularity. It may be noted here 
that if the series solution corresponding to different indicial 
roots are different, one can remove the term x m for negative 
and fractional m by the substitution y = x m u in (4*l) and carry 
out a similar analysis for the differential equation in u. It 
may also bo noted here that we need to transform the Chebyshev 
polynomial T*(x) valid for [0,1] to shifted Chebyshev polynomial 
T*(x) which are valid on (0,e] , e < 1 by the change of variable 


x = ^(x+l) 



Let us assume that for different indicial values the series 

Rl<,x) and ftjCx) are equal to R(x) so that the general solution 
is of the form 


m i m o 

y(x) = (a 1 x 1 + a 2 x 2 ) R(x) 


( 4*8 ) 


where m^ and rrt^ are the roots of the indicial equation* 

In order to approximate R(x) by an economized expansion 

N 

Pn(x), we assume that p^j(x) = S satisfies the differential 

r=0 


equation 


y"( x) +F^(x)y / (x)+F 2 (x)y(x) =t T*(x/e) (4*9) 

where T N (x/e) are the shifted Chebyshev polynomials in the 
interval 0 < x < s f and we choose T so that PjjCo) = a, where N 
and e are arbitrary constants* Now by substituting Pj^Cx) for 
y(x) and comparing the like powers of x on both sides in (4*9), 
we can find the coefficients b r and write the Eqn » (4*8) as 

m 1 m- 

y(x) = (a^ x + oc 2 x ) p N (x), 0 < x < e. (4*10) 


We now reduce the problem (4«l) ~ (4*3) to a regular boundary 
value problem, by finding a new boundary condition at x = e* 
To do this wo have from Eqn»(4*10), 


where 



y(x) ~ a 1 p(x) + a 2 q(x) 

ra 2 

p N (x) = p(x) and x P N (x) = q(x)* 


(4*11) ‘ 


Eqn * (4*ll) at x = e can be written as, 

y(e) = a 1 p(e) + a 2 q(e) (4.12) 

We also have from Eqn. (4.1l) 

y'(s) = p' (e) + a ? _ q'(e) (4.13) 

where primes denotes the differentiation-* Solving (4*12) and 
(4*13) for a 1 and ot^.* we get 

y(e) q'(e) - y'(e) q (e) 

1 p(e) q / (e)-p 7 (e) q(e) 

and 

y 7 (e ) p(s ) - y(e ) p 7 (s ) 

’ ' ‘ “ -* ■ 

p(e) q'(e) - p 7 (s) q(e) 

Since y(o) = a, we have from Eqn. (4*ll), 

p(o) 4 - a 2 q(o) = a 

using Eqns* (4*14), (4*15) and (4*16) we have 

y(e )q 7 (e )-y r (e )q(e ) y 7 (e)p(e) - y(e)p 7 (e) 

— • - p(o) + ■ ~ q(o) = a 

p(e)q / (e) - p 7 (e)q(s) p(e)q'(e) - p 7 (s)q(e) 

(4*17) 

or 

[ q 7 (e )p(o)~p 7 (e )q(0)] y(s )+ [p(e )q(o)~q(e )p(o) ] v'te) = ah(e) 

(4.18) 


(4*14) 


(4.15) 


(4-16) 


where 


h(e ) = p(e) q'(e) - p 7 (e)q(e) 



The Eqn* (4-18) can be 


conveniently written 


with 


and 


y(e) + Bj y* (e ) = C ± 

\ = q / (s)p(0) - p'(e) q(o) 
= p(e) q(o) - q(e) p(o) 

C t = ah(e ) « 


(4.19) 


(4.20) 


Eqn. (4*19) gives the new boundary condition at x = e • The 
reduced boundary value problem over [s/l] is given by Eqn*(4»l) 
subject to boundary conditions (4»19) and (4-3). 


In the case of non-homogeneous equation 

y' / (x)+F^(x)y / (x)+P 2 (x)y(x) = Fg(x) (4.2l) 

N 

the above procedure can be applied by making Pvr(x) = S b x r 

r=0 r 

satisfy the equation 


y # ( x)+P 1 (x)y / ( x)+P 2 (x)y( x) = \ T^(x/s )+F^(x) * (4.22) 

and obtain the coefficients b r , by comparing the coefficients 
on both sides of Eqn. (4*22)* 


4 » £ Invariant Imbedding : 

We employ the method of invariant imbedding for 
solving the regular boundary value problem given by Eqns. (4*1), 
(4.19) and (4.3). By reducing the boundary value problem 
into a set of initial value problems, we rewrite the Eqn. (4*1) 
as a systems of first order equation of the form. 
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u-'Cx) = v(x) (4-23) 

-v-'(x) = F^( x)v(x)+F 2 (x)u(x) 

The corresponding boundary conditions (4*19) and (4*3) 
as, 

u(e) + v(e) = 
u(l) =3 

Now let us consider a more general first order system given by 

u'(x) = a(x)u(x) +b(x)v(x) + e + (x) ( 4 * 27 ) 

x e [e r i] 

”'V / ( x) = c(x)u(x)+d(x)v(x) + e (x) (4*28) 

with all functions a(x), b(x), c(x), d(x), e + (x) and e (x) to 
be continuous in [e,l] - For the general system let us write 

u ( x) = S 1 (x) v(x) + S 2 (x)u(e)+S^(x) (4*29) 

x e [e,i] 

v(e ) = Q^(x) v( x) + Q 2 (x)u(e )+Q^(x) (4*30) 

It can be easily verified that the coefficients S ± (x) / Q ± (x), 
i = 1 , 2,3 satisfy the following initial value problems 

S^(x) = b(x)+[a(x)+d(x)] S^(x)+c( x)S^( x), S^(s) = 0 (4*31) 

S^(x) = [a(x) + c(x)S 1 (x)] S 2 (x) , S 2^ S ^ ~ 1 (4*32) 

S^(x) = [ a(x)+S^(x)c(x) ] S^( x) +S^(x)e(x) +e + (x), S 3 (e) = 0 

(4.33) 
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Q £(x) = [d(x)+c(x)S 1 (x) ] Q 1 (x) 

Q 1 (e) = l 

(4-34) 

Q^(x) = [c(x)Q(x)] S 2 ( x ) 

Q 2 (e) = o 

(4-35) 

Q^(x) = Q^(x)c(x) S^(x) 

o 

II 

CO 

a 

( 4 * 36 ) 


We evaluate (4*29) and (4*30 ) at x = 1* and solve (4*29), (4*30) 
and (4«25) as a system of three equations for the three unknowns 
u(c) wf v(e ) and v(l) * The solution u(x) for any x s is 

then derived by first finding v(x) from (4-30) and using it in 
(4*29) to get u(x). 

4*4 Computational Algorithm : 

To compute the solution* the computation is as follows : 

Step 1 ; Choose the values of N and e and determine the shifted 

Chabyshev polynomial Tj^jy'e), over the interval (0,e], 
s is near to singularity* Choose t such that p^(o) = a 
and find the coefficients b r in the polynomial p^Cx) by 
solving an algebraic system resulting from making p^C x) 
satisfy the Eqn* (4*9)* 

Step 2 : Find the boundary condition (4*19) by computing the 
values of A^, and C^ * 

Step 3 : Integrate the initial value problems given by Eqns* (4*3l) 
(4*36), using efficient initial value routines from 
x = s to x = 1 to obtain the profiles S^' s and Q^'s and 
store them* 

Step 4 : Evaluate the Eqns* (4*29) — (4*30) at x = 1* Using the 
values of S ± ( l) and Q^l), (i = 1/2,3) from Step 3 and 
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u( l) from (4 .26) . 

Stop 5 i Combine the equations from Step 4 with Eqn* (4*19) 

(from Step 2) and solve these three equations for the 
unknowns u(e), v(e) and v(l). 

Step 6 : Compute the other values of v(x), u(x) for x £ 

by using the values of u(s)., v(e) from Step 5, and the 
stored values of S ± and Q ± (for i = 1,2/3) from Step ?* 

4*5 Numerical Results and Discussion : 


In this section the numerical results of two model 
problems solved are given to demonstrate efficiency of the 

method* An computations were carried out on DEC-1090 Computer 
system- 

Puoblem 4*1 : We consider 


2x( 1+x) y" + ( 1+5 x) y 1 +y = 0 
y(o) = l 

y(l*5) = 1 * 


This problem has earlier been studied by Cohen and Jones [26] 
by using finite difference and deferred correction approach* 

The analytical solution for the above problem is given by 


y( x) 


l-Wi *5x lj (f 

— ( l+x')"' • 

The polynomials for N = 5 and e =0*5 and 0*1 are given 


by y 
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P N (x) = 1 #0-0 *9999 1287x+0 •9846067x^-0 *891666 1x^+0 ♦6l28445x 5 

-0.216 298 lx 6 

and 

p N (x) = 1.0~0.9999991x+0,9992441x 2 -'0.9968280x 3 -+0.9605363x 4 

- 0.7O2366lx 5 

respectively. 

Problem 4.2 : We study another problem 

4x(l+x)y w + (?+llx)y / +y = 0 
jiO) = 1 
y(l«5) = l. 

The analytical solution for the above problem is given by 
y ( x) = l+(l .5) 9 / 8 x 1//4 /1+x» 

For our numerical experiments, we take N =8, e =0.5 
for which the polynomial p^Cx) is given by 

Pjj(x) = 1*0-0 .9999985x+O.9999175x 2 ~0.9984735x 3 +O.9857662x 4 

-0 .9 227079 x 5 +0.7 3749 32x 6 -0 *4178441x 7 »0 . 118024 3x 8 . 

The numerical results for Problem 4*1 for N = 5 and 
e = 0-5 and 0.1 are given in Tables 4.1 and 4.2 respectively. 
The computed solutions at various values of x compare very well 
with the analytical solution* This example has earlier been 
considered by Cohen and Jones f26j, who have solved it using 
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finite difference deferred correction technique* They have 
used economized series expansion in the interval [o,lJ and 
obtained finite difference solution on the remaining part of 
the interval - By doing so they neglected the effect of 
singularity on the solution in the immediate neighbourhood of 
the singular point since the difference solution computed is 
far away from the singularity* Their solutions are also 
presented in Tables 4-1 and 4-2 and it is observed that the 
solutions in the viscinity of the singular point computed by 
invariant imbedding method are quite comparable to those of 
Cohen's solution obtained from economized expansion there* It 
has also been observed that for £ < 1, a comparative smaller 
value of N is required than for the case when 6=1, for achieving 
the desired accuracy of the solution- Its advantage lies in 
solving a comparatively smaller set of algebric equations to 
compute the coefficients of the polynomial p^(x)- This saves 
on both computational time and memory requirements and also 
minimizes loss of accuracy due to rounding errors* The computed 
solution which are presented in Table 4-2 are seen to be 
matching with the analytical solution upto atleast five places 
of decimal- The numerical results for Problem 4-2 for N = 8 and 
s =0-5 are given in Table 4*3 and agree with the exact solution 
upto five places of decimal- 



TABLE _ 4,1 



Numerical 

Results for the 

Problem 

4-1 


(N 

= 5 and s = 0 *l) 





*"'** """ ' X * ' f * r - 

*“ - m: - ■* •• 


X 

Imbedding 

Solution 

Cohen' s 
Solution 


Exact 

Solution 

0 .1 

1*261176 

1*261183 


1*261180 

0.2 

1.239774 

1*28977 2 


1 *289769 

0 .3 

1.28525 2 

1*236251 


1.285247 

0 *4 

1*267575 

1*267574 


1*267569 

0 .5 

1*244023 

1*244022 


1 *244017 

0 .6 

1*2179 30 

1*2179 31 


1-2179 27 

0.7 

1*191003 

1*191000 


1.190997 

0 .8 

1 *164140 

1*164139 


1.164136 

0 *9 

1*137844 

1*137842 


1-1378 39 

1 .0 

1*112362 

1 *112 374 


1*112372 

1 .1 

1 *037370 

1*037870 


1 .087368 

1 *2 

1*06 4 384 

1*064303 


1 *064382 

1 *3 

1 *041923 

1 *041924 


1.0419 23 

1 *4 

1*020474 

1.020474 


1*020474 

1 ,5 

1 *000000 

1 *000000 


1 *000000 
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TABLE 4.2 

Numerical Results for the Problem 4-1 
(N =5, and e =0.5) 


X 

Imbedding 

Solution 

Cohen 7 ' s 
Solution [26] 

Exact 

Solution 

» .to. hr ' 

n.-BTW -■»:*. » - *■ ktx w—' »- l. - 


■ — ■ - - ■ — - 

0-5 

1*244025 

1*244022 

1 .244017 

0 *6 

1-2179 33 

1*217931 

1.2179 27 

0*7 

1-191010 

1 *191000 

1.190997 

0*8 

1*164138 

1.164139 

1 .164136 

0 .9 

1 -137846 

1*137342 

1 -1378 39 

1*0 

1-112376 

1 *112374 

1*112372 

1 *1 

1*087372 

1.087870 

1 .037368 

1 *2 

1.06 4 384 

1*06 4 38 3 

1 -064382 

1 .3 

1*0419 26 

1 *0419 24 

1.0419 23 

1 *4 

1 *020474 

1*020474 

1 .020474 

1 .5 

1 .000000 

1 *000000 

1 .000000 


, _ 

. ^ ♦ *. .. ,* „ , , 'm . » 

. . • • • * • 



TABLE 4*3 


Numerical Results for the Problem 4*2 
(N = 8, and e =0.5) 


X 

Imbedding Solution 

Exact Solution 

0.5 

1.4265032 

1.426502 

0 .6 

1.3705673 

1.370566 

0.7 

1.3175152 

1-317515 

0 .8 

1.2677009 

1.267701 

0 *9 

1.2211408 

1.221141 

1 .0 

1.1777007 

1.177701 

1 »1 

1.13718 38 

1-137184 

1 *2 

1.099 3700 

1.099 369 

1.3 

1.06 40368 

1.064037 

1 «4 

1.0309788 

1.030979 

1 .5 

1.0000000 

1 .000000 



CHAPTER 5 

MODIFIED FOURTH ORDER FINITE DIFFERENCE 
METHOD FOR LINEAR SINGULAR BOUNDARY VALUE 

PROBLEMS 

5*1 Introduction: Most of our attempts towards the development of 
numerical techniques for linear singular boundary value problems 
mentioned in the previous chapters are based on the removal of 
singularity by making use of ejqpansion procedures in the vicinity 
of the singularity- However the method of finite differences 
is another practical method suited to the numerical solution of 
singular two point boundary value problems * However for some 
problems, it may sometimes be difficult or even not possible to 
obtain the series solution in the neighbourhood of the singularity <> 
In such cases, one may need a direct method to solve singular 
boundary value problems* The method of finite differences, which 
is quite popular for solving boundary value problems, can be one 
such direct approach, to tackle these problems* 

Finite difference methods fob singular problems have been 
studied by several authors namely Jamet [ 50] , Gustaff sson [ 47] , 
Chawla [2l], Deodel and Reddien [37], Russell and Shampine [83], 
deHoog and Weiss [29] , and Brabston and Keller [16] • Chawla [2l] 
has discussed the construction of three point finite difference 
approximations and their convergence under appropriate conditions 
for a class of singular two point boundary value problems* This 
has been done by establishing a certain identity based on a 
general mesh, from which various methods are derived* Also 
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Russell and Shampine [8 3] have used the traditional second 
order finite difference, patch bases and collocation methods 
for singular boundary value problems* Recently Deodel and 
Reddien [ 37 ] have constructed some finite difference schemes 

these problems and showed how these methods can be viewed as 
a special type of collocation scheme* 

The object of this chapter is to present a modified 
fourth order finite difference method to solve a certain class 
of linear singular boundary value problem* 'The original 
differential equation is modified at the singular point- The 
main feature of the modified difference scheme is that it leads 
to the tridiagonal system of equations which has been treated 
by discrete invariant imbedding method* The Numerical experiments 
for the model problems have been given to illustrate the method. 

5*2 Difference Method : 

We consider a certain class of singular boundary value 
problems given by 

— -•■X +. — OX - q(x)y(x) = f(x), 0 < x < 1- (5*1) 

dx x ax 

with K > 1 and subject to boundary conditions 

y' ( 0 ) = 0 ( 5 *2 ) 

y(l ) = p . (5 * 3) 

As Jamet [50] has shown that for the Eqn* (5*1), the derivative 
boundary condition is irrposed due to nature of the physical 
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situation of the problem# Due to singularity at x = 0, we need 
modification of the problem near the singular point* We present 
a modified fourth order finite difference methods which leads to 
a tridiagonal system* 


We divide the interval 0 =x q <x 1 <x 2 = 1 intc 

n equal parts with uniform mesh size x^— x^^ = h, where h is. the 
mesh size* By using a fourth order difference method and 
employing central difference formulae for the first and second 
order derivatives [ 42 ] , we may write Eqn* (5*1) as. 


{<5‘ 


~ U 6 ' 


JL 

90 


6 - 


1 

560 


8 




iY ± 


+ hxT t 116 ~ k ^ + + 

~ q^y.^ = ~ 1/ 2, * * »/n— 1 

where y± = yCx^)* 

Rewriting Eqn* (5*4) in a compact form. 


(5*4) 


-=7 i<5 2 - ^7 6 4 } Yi + - h "x“ ~ h ^ 3 }y i -q i y i = f i +E ± ^5*5) 

h * . i 

where E i = ^ +•»*} Y± - £x. *55 V6 5 +* } Y ± (5*6) 

h . i 

and E ± = 0(h 4 ) (5*7) 

Eqn* (5 *5 ) is a fourth order finite difference method which is 

3 

not tridiagonal due to the presence of the differences y ± 

and in the left hand side of the equation* To make 
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Eqn. (5.5) a tridiagonal system, we ^ 


of U6 y* and 6% upto orfer f±ve ^ 


the tridiagonal estimate 


six respectively. 


he singular point x = 0, we modify Eqn.(5.1) as. 
(1+K)y"( x ) _ g(x) y (x> = f(x)< x = ^ (s<8 

FOr Bln ' (5 ’ 9) ' we a PPly the above mentis _ .. 


x = x _ 


( 5 *8 ) 


difference formulae at x = x ( v _ n \ ^ ^ 

^3 ’ x — 0 ) and obtain 


mentioned fourth order finite 


iJ - rx2 l 4 


^ - f6 — 1 a 4, 

r ■ n ! y o - OoYo - f 0 +H 0 


(5.9) 


where H = f L „ JL 

° 1 p * 90 


k . JL 5 6 1 6 ; 

° + —n* » 


* 560 '*** Y 0 


Differentiating Bgn. (5*) twice with respect to x yields. 


(l+K)y iv ( x )^ q *( x ) y(x )^ 2g / (x)y/(x) ^, (x) = ± *( x ) f x = x 


Y ~ Tl+kT fq(x)3r^x)+2<i'(x)y / (x)+q*(x)y(x)+f' r (x)} / x = x 


y (x Q ) = ^ ^' +K J { g( 3 ^ ) y*( x 0 ) } + J- q , (x o )y / (x o ) 


q"( Xq ) f "(x) 

+ ri+io Ylx o> + ri+Kr 


(5-10) 


4 ^ 3 

6 y(x o ) = (i+k 7 < 3 (3 S ) 6 2 y(^) + fy +l0 g /(x o } M6y(x o ) 


hV^^) h 4 f«( x ) 

+ ilffl + • U+K") - .;+ G ± 


where Gy = O(h^) 


(5.11) 
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Substituting Egn. ( 54l) ^ (s>9) 


we get 


s 2 y - i.l+K) ,h 2 
v 2 y o . 2 " ir'i 


12h 2 " { ^«) q(x o } 5 y(x o ) + T ^ kJ . q'(x ) 


M6y(x Q ) 


h 4 q"(x ) h 4 f- a ’( v ) 

+ -Ti + S y(x D ) + _. >v 


By neglecting the truncation 


we get 


1+k; ^Vo^i = f i (5.12) 
where A . - 0(h 6 ) 

error term and after some algebra. 


r-t 

L 1 1 2CU 


h 3 


'*£> . r 5h \ h 4 q£ 


2 Ci+K J "* i'2(i+k)1|*^6'(i+K7 + 12 T 1 +KJ + 2 J y ( 


i_2 , 3 

h q o . h X 


o 

+ fi -o “ 4 q h f 4 

L ~ 1 2 ( i+k 7 + i 2 (i +K ) ]y_i = U-HO + h f " 

(5.13) 

where = q(^), = q'(^), f Q = fC^) and f' = 

Using Eqn. (5.2) and (5.13), Ekqn. (5.13) can be written in a 
compact form. 


where. 


(l» + N )y ~ M y = p 

0 O -*1 o^o 0 

(5.14) 

•u 3 , 2 

L V ^ 

O " 12( 1+k7 “* l2n+k) + 1 

(5.15) 

5h^q h 4 q* 

M o = 6 ( 1+K J + lYC l+KO + 2 

( 5 .16 ) 
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and 


N 


h^rr' v>2 

_ ■%> h % 

12(1 +K 7 ~ i 2( i+Kj + 1 



~ Ti+k) 


hf 
+ 12 



(5.17) 

( 5 .18 ) 


Differentiating Eqn»(5»l) with respect to x, we get the tridiagonal 
estimate of P*6 3 y ± (i = 1, 2, . . .,n~l) as 


T tn 


K 

x 


~K 


Y± ” T *± ~ ^i^i ~ = f i / (i = 1,2, 


or 


c 5 hIC 2 Kh - ,,, , 3 . 2 3 , 

y ’ 6 Y i + xj 6 Y i ~ h g5.Yi ~ ft q ± M6Y i = h 


(i = 1,2/ • • n-1 ) 


or JLlfi^y . = - ~ -2 ,2/K n 

i x i 6 Yf + h (-j + g ± ) H6y i 

x i 

+ h 3 q^y ± +h 3 f£ + 0(h 5 ) (5.19) 

Similarly differentiating Eqn*(5»l) twice with respect to x, 

4 

we get the tridigonal estimate of 6 y^_ as 


s^y^ = h 2 (^?r + + q-j_) s 2 y± 

x i x i 

+ h 3 {2q , + ^i.2K.2K? )M6yi . 


X j 


Kq, 


+ h 4 (q » + |- ql - £± ) 


+ 0(h) 


(5.20) 
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Using Eqn * (5*19) and (5.20) into (* k \ * ■ 

; into iqn* (5*5) and rearranging the 

terms we get a tridiagonal system 


where 


and 


L i 

y itl 

- M ± y 

i + N i y 

i -1 

p i' 

i = 

= 1,2 

/ *••/ 

n— 1 ( 5 . 21 ) 

L i 

= 1 ■ 

. Kh 
25- 

+ — ( k2 
+ T5 """2 

x i 

K 

~ T 

x i 

~ q ± ) 4 


,2K 

^ 3 “ 

x i 

Kq. 

24 - 

i . 










(5*22) 

M ± 

h 4 

= rs 


K ■ , 

^ q i~ 

K N 

T q i } + 
x i 

h 2 

6 

(5q ± 

+ K 2 

+ T 

x i 

-^) + 2 
x i 










(5 .23) 

N i 

= i~ 

hK- , 

+ 

h 2 f K 2 . 
x i 

K 

T ~ 

x -f 


‘ 24 

, 2 K 

- 24 

x i 


’± = h ‘ £ i + (f ± + i: f i + t f ±) 


K- 


(5*24) 


(5.25) 


Eqns. (5*14) and (5«2l) contribute a difference system of n 
equations with n unknowns. 

5 *3 Computational Method : 

To solve the tridiagonal system we make use of the 
method of invariant imbedding [ 7] * We seek a difference relation 
of the form, 

Yi = Wiy i+1 + Tj i = 0,1, 2 , . *.,n— 1 (5.26) 
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where W ± and T ± correspond to w(x A ) and T( Xl > and are to be 
determined * 

Using Eqn. (5*26) in (5,2l), we have 


or 


or 


L i y i +1 - Vi + B i<Wi_iyi + - p i 

L i Y i+l + < "i K i-l" K i ' 1 j'i = p r N i T i_i* 

Y i [ H i w i-r M i] - - L i Y i + i + C p i- H i T i-i] 

L i £ p i - Vi-il 

y i - TM^w-r y i + i + 


(5*27) 


Comparing Eqn* (5*27) with Eqn* ( 5 *26 ) # we see tha t the recurrence 
relations for and T^ for i = l/2*»**,n-l are obtained as 


W. = 


(M ± - N i W i _ 1 ) 


and 


_ p i.~ ¥w_ 

1 ’ < N i w i-l- M i ) 


( 5 .28 ) 


(5.29) 


To solve the recurrence relations for and T^, we need to 
know the initial conditions for W Q and T q * These initial values 
can be easily obtained as. 


W_ 


(L 0 + N 0 ) 
M. ’ 


(5.30) 
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rn _ U 

o'- M~ (5*31) 

o 

where L q/ H q/ N q and P Q are given by Eqns. (5.15) - (5.18). Using 
the initial values, we compute VI ^ and Th for i = 1,2, »*.,n-l / 
from Eqns. (5.28) and (5.29) in the forward sweep manner and 
then obtain the solutions y^ in the backward process using the 
Eqn • (5.26). 

5-4 Numerical Results and Discussion : 

In this section the numerical results of the model 
problems solved are presented 

Problem 5 *1 : We solve 

Y" + | y' - 4y = —2 
y* ( o) = 0 

y(l) = 5.5 

This problem has earlier been solved by Russell and Shampine [83], 
which has the exact solution 

/ \ 15 sinh 2 x 

y(x) = 7 + - x stoh ‘2 

Problem 5.2 : We solve another model, a non— homogeneous equation 

~y* „ | y' + ( l-x 2 )y = 7 - 2 x 2 +x 4 
y* ( 0 ) — 0 
y(l ) = 0. 
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This is the model problem used by K. Eriksson and V. Thomee [39] 
and has the exact solution y(x) = 1-x 2 . 

The computational results for the Problem 5*1 for 
different mesh sizes are given in Table 5*1* It is evident from 
the results that the modified finite difference method works 
well and the solutions are quite comparable with the exact 
solution* The results indicate that near the singularity one 
may not need a very small mesh size due to the modified scheme* 
The maximum error between the approximate solution and the 
computed solution at the mesh points for different values of h 
are given in Table 5*2* Also comparison of the maximum errors 
with regard to other methods is presented in Table 5*2* It is 
observed that the fourth order finite difference schone gives 
better accuracy than other methods* 

The results for the Problem 5*2 are given in Table 5*3* 
The results indicate the efficiency of the method for the non~ 
homogeneous problem also* 
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CHAPTER VI 


DIFFERENCE METHOD FOR NONLINEAR SINGULAR PROBLEMS 

6«1 Introduction: The object of this chapter is to extend the 
finite difference method described in Chapter V to non-linear 
singular boundary value problems. We shall consider a class 
of non— 1 in ear problems of the following form : 

y /# (x) + y*( x ) - q(x)y(x) = f(x,y) (6.1) 

subject to the boundary conditions 

Y'(O) =0 (6.2) 

Yd) = 3 (6.3) 

where m > 1 and |3 is given constant. We assume that (i) f(x,y) 

is continuous ( ii) ~ exists and is continuous and — >0. 

o y ay - 

The application of finite difference methods to these problems 

has been discussed by some authors namely Jamet [ 50], Deodel 

and Reddien [37] * A three point finite difference scheme has 

been proposed by Russell and Shampine [8 3] with Newtons 

iteration procedure for solving non-linear problems* 

In this chapter we have presented a method based on 
finite difference technique for non-linear singular boundary 
value problems * The quasilinearization technique, originally 
developed by Bellman and Kalaba [12], has been used to reduce 
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the given non-linear problem to a sequence of linear problems# 
The fourth order finite difference method, similar to the one . 
described in the previous chapter, has been applied to solve 
the linear problems# Some model problems have been solved to 
demonstrate the efficiency of the method. The comparison of 
numerical results is also made with other methods# 

6 #2 Quasilinearization In the Quasilinearization technique* 
instead of being solved directly, the. non-linear differential 
equation, is solved recurssively by a sequence of linear 
differential equation?# The main advantage of this method is 
that if the procedure converges* it converges quadrat ieally 
to the solution of the original problem# Quadratic convergence 
means that the error in the (K+l)th iteration is proportional 
to the square of the error in Kth iteration* The linear 
equation is obtained by using the first and second term in the 
Taylor's series expansion of the original non-linear differential 
equation* 

Consider Eqn* (6.1) with boundary conditions (6.2), 
(6.3)# We choose a reasonable initial approximation for the 
function y(x) in f(x,y), call it as y^ (x) and expand f around 
the function y°(x), we have 


f(x,y (l) > = f(x,y (o) > + (y ( 1 , -y (o> > HI, (o> ) + " 


(6.4) 


Or, in general we can write for k - 0,1,2, ••• 
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) = f(x,y <k) ) + ( y Oc+l> (k)) St. 


3 y: (x,y (k> 


Eqn* (6*1) can be written 


as 


(k+l) + m , ( k+i) ^ 


q(x)y (k+l) = f(x,y (k) ) 


+ (y (k+ l> - y< k >) U 


3 y (x,y (k) ) 

k — 0,1,2^*#* 

subject to the boundary conditions 
y' (k+l) (o) = o 
y (k+l) (l) = 3 

Let us re-write the Eqn* ( 6 * 6 ) in the following form 


y* (k+l) + 5 Y' (k+l) - r (k) (x)y (k+l) = F (k) (x) 

where 

(k) , n / w .3f 
r (x) = q(x) + 3 -y ( k ) 

J ( x, y J 


and 


P (k) (x) = f(x,y^ k ^) 

k 0/1,2/ * * * • 


(k) 


afi 

3yi 


, (k) v 

(x,y ) 


6*3 Difference Method : - In order to solve the system 
linear singular boundary value problem given by Eqn* (6 
subject to Eqna * (6*7) — (6*8), we use the fourth order 


) 

(6 *5) 

(6 *6 ) 

(6.7) 
(6 * 8 ) 

(6 #9) 

( 6 * 10 ) 

(6 * 11 ) 

of 

►9) 

finite 
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difference method described in the previous chapter. Due to 
the singularity at x - 0 we need to modify the problem at x = 0. 
As usual, we divide the interval [ 0 ,l] mto N equal sub- 
interval with fixed step length h such that o =x Q < x, < ... < 3^=1 
By employing central difference formulae for the first and 
second order derivatives and taking upto fourth order terms we 
get the difference equation for (6.9) as 


“V “ 1 2 Yi k+1) + [M6- £ U6 3 ] 

- r^ y^ k+1 ^ = + E ± 

where = r^ k \x i ) and F^ = F^ k ^(x ± ) 

and 


<k + n 


and y^ k+1 ^ 


(6 # 12 ) 


= y (k+1) ( Xi : 


Ej = 


1 


6 6 +*•*} 


y' k+l) 


m . -j 
hx, 1 30 


46 5 +.*.} y k+1 (6*l3) 


with i = 1,2, •*#,N-1# 
and k = 0,1,2, *##« * 

As discussed in Chapter V, at the singular point 
x = Xq, the modified differential equation can be written as 
follows 

( 1-Hn) y* (k+l) (x) - r (k) (x) y (k+l) (x) = F (k) (x), x = x Q , 

k = 0/1*2/*** • ( 6 *14 } 


Re-writing Eqn* (6 #14) we have 


A 
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" ( k+1 } ( x) - ^ v (k+l ) ( ^ _ F ( k) ( x) 


TfSnT y 


(v - ) _ £ vx. 

lx) ~ Ti-taT 


(6.15) 


We will now employ the fourth order finite difference method at 
x — 0 for tnis modified equation (6*15). Then we get/ a 


difference equation at x = x Q (= 0) as 


L r 6 2 „ X. fi 4 l v (k+l) r (k) (x) (k+l) *o"' . „ 

2 L 12 5 J y o ~ Ci-mT ~ y o = Ti«T) + E o 

(6,16) 

6 


,(k) 


where 


E — — -.--L f ™ 4. 

o ~ h 7 L 9 0 + 


•] yi* +1) 


To make the equations (6*12) and (6*16) as tridiagonal 

3 ( ) 

system, we first make the tridiagonal estimates for H 6 y^ , 
6 4 y^ k+ ^^ and 6 ^y^ -1 "'*' ^ as follows* Differentiating Eqn« (6*15) 
twice we get at x = x D , 


T ru 


<k+l) (x) - y (k+l) (x> - y' Ck+l) (x) 


F /(k) (x) 


X 


*o 


r iv(k + l) U ) _ y(k+l) (x) _ y" (k+l) U> 

_ r^UX y -(k + D (x) . 4^. y» '*«’(*) 


J// (k)( x ) 


x = x,_ 


(6-17) 



101 


The tridiagonal estimate of 6 4 y^ k+1> (* = 0,1,2,...) for 
Eqn. (6 #17) yields. 


4 ( k+l) 


6 Y, 


. r (k) , v ~*(k) 

h ( Tf«T )4 V 0 ]H1) + h 3 ^-) n 6Y ^«-> 


+ h4( Tv™r ) yo K+1 ' - » 4 (^. r ) + o(h 6 >. 


r? (k) /F'*(k) 


(k+l) , .4, o 


k = 0,1,2, *»• ( 6 *18 ) 

Substituting Eqn« (6*18) in (6*16) and neglecting the error terms.; 
after some manipulations we get, 

« (k) -s _/(k) 

n h ° _ h r o .(k+l) 

L 1 rz Ti+mT rs n+mT ] y i 


2 r (k) 4 v/» (k) . 

To A 5h f o v . h* ,*0 \i __( k+l ) 

- L 2 + "6 1 ( Tl+m) } + 12 (T+mT’ -* y o 


o Jk) - *(k) 

r, h 2 r o h 3 *Q (k+l) 

+ Ll ~ 12 n+ml + 12 (T+m) J y ~l 

h 2 p^ k ^ .4 F M ^ 

= Xl+mr + n Ti'+m)' k = 0 , 1 , 2 , 


(6-19) 


By using Eqns. (6*7) and (6*19), we write Eqn* (6*19) in a 
compact form as 


(k+l) _ n v (k+l) _ £, 


( A o + c 0 ) - B o y o 


where. 


(6 . 20 ) 
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and 


A_ = 


B_ «e 


C_ = 


D = 


2 r (k) - (k) 

1 _.L. lo. _ h J ^ 

12 Tl+m) ~ 12 Xuaj 

(6.21) 

, h 2 r (k) .4 r" (k) 

6 Tl+m) + 12 TT-wi J~ 

(6.22) 

h 2 r< k) h 3 r^ (k) 

n I l+rn*J + 12 Tl+mJ 

(6 -2 3 ) 

h 2 p0<) ^4 (k) 

‘tvhSt + T 5 ifsrr 

( 6 - 24 ) 


for k — 0,1,, 2/ * 

Again differentiating Eqn* (6*9) twice with respect to x, 
and neglecting the error tenns, we get the tridiagonal estimates 
of and 6 4 y.j_ (i - 1,2, «* >,N-1 ) as. 


^6 3 y ( i k+l) 


h 2 (k+l) .2, H ( k) s ... (k+l) 

= ” — 5 y i - h l ~ r ± ' ^y± 


and 


+ h 3 r' 1 ( * ) y ( 1 * +l) + h 3 ]^ (k) + oCh 5 ) 

k = 0/ 1 /2/ * * * 

r 4 , 2,m 2 . m . Jk)\ ,2 (k+l) 

6 y^ = h \ — 7 £ + — ^ 4- r\ / 5 

x i x i 

+ h 3 (2r' (k) + ~~ 4 k) ~ ^ ~ x" } M6y i 1C+l) 


+ h 4 (rf' (k) + f- ri (k) 

1 x 


x . ' i 

1 (k) 

mr i ) _<k+l> 
- * ■) y ± 


+ h 4 (F" (k) - £- Bl (k) 4 k>) + ° C1,6) 

1 X i X i 

k = 0,1,2, 


(6-25) 


(6-26) 


* • m ; . # 
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Substituting Eqns. (6 .25) - (6.26) in (6.12) and rearranging 

the terms we get a tridiagonal system (by neglecting the error 
term E^) as, 

A i y i+i l) - B iy ( i k+1) + = Dl , i = i, 2,...,N-a 


k =0,1,2, — . (6.27) 


where 


A - 1 . k m , ^ t XU ' 

i- 1+ T37 + T2 ( T~~I 

i x ± x 1 


- r 


(k) 


, h .. 
+ 72 ^ 


(k) 

(2m « ,(k) i -v 

1 X i 


X . 


( 6 . 28 ) 


= 2 4* ( 5 it 4* 


h' 


2 

m 

m 

T " 

' “7 

x i 

x i 

(r '! 

(k) 


. (k) (k) 

mr\ mr. 

+ --£r- - i--> 


( 6 .29 ) 


and 


, h m 
c i = 1 ~ 7 


h ? ‘ ,m 2 m (kK 

x, + n ( ;t ~ ^ ~ r i 3 

k i x i 


„ (k) 

_ £ (2S - ar A *> - 3 L .■-) 

13 ( j 2r i Xl 


^ ^ ^ r i !k> 4 l,) 


(6 .30) 


(6.31) 


where r (k) (x) and F (k) (x) are defined by Eqns. (6,10) and (6.1l) 


respectively. 
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We observe that Eqns. (6.20) an d (6.27) produce a 

tri-diagonal system of n equations with n unknowns, which can 

toe solved by the methods discussed earlier* The approximate 

solution values, at each stage of the iteration (k = 0,1,2, .**) 

are computed till the convergence criteria iv^ + '^ — < EPS, 

2 i 2 i 7 

where EPS, is given value, is satisfied. 

6 .4 Numerical Results and Discussion : 

In this section, some physical model problems solved are 
presented* All computations were done on DEC~1090 computer 
system with single precision arithmetic. 

Problem 6*1 ; We consider the problem arising in Astronomy, 
the equilibrium of isothermal gas spheres can be described by 

y"(x) + ~ y'(x) + y 5 =0 

with boundary conditions 

y'(o) =0 

y( 1 ) = '/IT?- 

*2 

which has a exact solution y(x) = l/'fl + ^x . This problem Is 
stated in Russell and Shampine [83] , Hoog and Weiss [29] and 
Rentrop [ 53 ] . 

Problem 6*2 : We solve 

y"(x) + -. y'(x) + e y(x) =0 

A 
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with boundary conditions 

y'(o) = o 


y(l) = o 


This problem has earlier been discussed by Russell and 

Shampine [83] • The exact solutions are y(x) =2 log(- B -+ 1 . ), 

Bx 2 +1 


where B = 3 + 2^2* 


Problem 6*3 : The equation for one dimensional heat conduction 
with non-linear heat generation is given by 

Y"(x) = ~ J y'(x) -Xe y(x) 

A 

with boundary conditions 


y' (o) = y(l ) = 0 


where X t heat generation constant, O < X < 0*8 

y : temperature distribution* 


This problem possesses a numerical singularity at x = 0. This 
problem has been considered earlier by NA [6 7] • 

Problem 6 *4 % As a last example, we solve 

? rp(l~y(x)) 

y^x) + ~y'(x) = f y(x) exp [ - - --- - - ] 

x 1 +(3(l~y(x)) 


subject to boundary conditions y^Co) - y(l) = Oi 


y % dd.m 0 nsi.onl 0 s s cone ont its. t don 

a : structure of catalyst particle, a = 2 spherical 


j 3 ,7, <p ; chemical constants* 
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This problem is studied by Rentrop [67] , for different values 
of J3, V and <p. 

The numerical results for problem 6.1 for different 
values of mesh size h are presented in Table 6*1. These solutions 
correspond to the iteration index k = 4 by taking = 0 as 
the initial approximation. The iterations were stopped by 
satisfying the absolute error criterion lyj_ k+1 ^-yjj*^ 1 < ICT 5 for 
it is apparant from the table, the results converge 
to the exact solution. The comparison of the computed solution 
of this program for h = 1/64 with other methods is given in 
Table 6.2* It can be observed that our solutions compare well 
with the solutions obtained by other methods. It can also be 
noticed that the solution near the singular point is not affected 
much due to the removal of the singularity at x = 0. 

The numerical results for Problem 6 *2 for different 
value of h are given in Table 6 ..3* These results correspond 
to the iteration index k = 3 and y Q = 0 as the initial 
approximation. The comparison of the computer solutions obtained 
by different methods is presented in Table 6.4* It can again 
be observed that the results agree well with the smaller 
solution, which is the exact solution with the value of B equal 
to 3-2 V2. The numerical results for problem 6.3 for X = 0*8 are 
given in Table 6.5* These solutions correspond to k = 4 and 

y Q = 0» 
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The numerical results for Problem 6.4 are given in 
Tables 6*6 and 6*7» Table 6 *6 presents computed solutions for 
(3 = y = (p - i and the computed solutions for ,8 = 0*05, Y = 20 
and (p = 6 are given in Table 6»7. The value of a in both 
these cases has been taken equal to 2. These solutions 
correspond to y Q = 0 and k = 4 satisfying the stopping criterion 


( k+1 ) (k), . , _-6 - , ■ . 

r i - Y. , I < 10 for all i. 


The solutions are seen to compare well with the solutions 
obtained by Rentrop [67] • 
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Table 6.1 

A " . r a — r 

Numerical Results for Problem 6.1 


X 

10 

20 ' 

40 

60 

0.0 

0.9976150 

0,999 3991 

0.9998491 

0.9999 315 

0.1 

0.9976565 

0.998159 2 

0.998 2917 

0.998 3156 

0-2 

0.99 27323 

0.9932253 

0.99 33547 

0.993 3780 

0 *3 

0 .9846899 

0.98516 31 

0.98 5 2869 

0*9853090 

0.4 

0.9737577 

0.9742001 

0.9743154 

' 0.9743359 

0.5 

0.9602 305 

0.9606 300 

0.9607338 

0.96075 21 

0.6 • 

0.9 444484 

0.9447922 

0.9443813 

0.9448968 

0.7 

0.9267758 

0 .9 270509 

0.9 271218 

0-9271342 

0.8 

0 .9075811 

0 .9077749 

0.90782 47 

0.9078 335 

0.9 

0.8 8 72205 

0 *8873219 

0.8 873479 

0.8873525 



Comparison Table for Problem 6,1 




,r - — - - * « 

" " - - * * - 


X 

Patch 

bases 

[83] 

Finite 

Difference 

[83] 

Our 

Method 

Exact 

Solution 



- ,r_ -r . w , 

— * “ " » ** ' 

- - ' — - — 

0 ,000 

0.9999 2 

1 .00000 

0.9999 372 

1 .OOOQOOO 

0 ,1 25 

0.99737 

0.99742 

0.9973843 

0 ,99 74060 

0 ,250 

0 .989 71 

0.98976 

0.98972 31 

0,9897433 

0,375 

0.9 7733 

0.97737 

0.9773373 

0.9773556 

0*500 

0.96075 

0.96078 

0,9607530 

0.9607689 

0,625 

0.94062 

0 ,94064 

0.9406 213 

0.9406342 

0 ,750 

0.91765 

0*91767 

0.9176537 

0.9176629 

0.875 

0.89257 

0,89256 

0.8925647 

0.8925696 


Tab! e 6 *3 

Numerical Results for Problem 6*2 


X 

10 

' * ‘Vr r , -.r . . 4 *• 

20 

40 

60 

0 .0 

0.3161375 

0.3165530 

0.3166751 

0 .3166751 

0*1 

0 .31 289 75 

0.3131708 

0 .3132412 

0.3132533 

0.2 

0 .3026895 

0.3029 309 

0.30299 34 

0.3030041 

0.3 

0 .2857503 

0.2859701 

0.2860272 

0.28 6 0369 

0 *4 

0.262 2603 

0 *2624607 

0.2625128 

0.2625217 

0 #5 

0.2 3245 34 

0.2 326 335 

0.2326803 

0*2326884 

0 .6 

0*1966152 

0.1967719 

0.1968125 

0*1968197 

0.7 

0.1550750 

0.1552033 

0.1552 365 

0*1552424 

0-8 

0.10819 6 7 

0.1082902 

0.108 3143 

0.108 318 7 

0 .9 

0.056 3698 

0.0564209 

0 *0564340 

0.0564364 
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Table 6 *4 

Comparison Table for Problem 6.2 
(h = ^ 



Finite 

Patch bases 

Our 

Smaller 

X 

Difference 

[83] [83 J 

Method 

Solution 

0 *000 

0.31762 

0.31672 

0.3166778 

0.31669 

0.125 

0 -31135 

0.31135 

0.3113296 

0.31134 

0.250 

0.295 37 

0.29537 

0.2953528 

0.29536 

0.375 

0 .26909 

0.26909 

0.2690043 

0.26901 

0.6 25 

0.18696 

0.18696 

0.1869479 

0.18695 

0*750 

0 .13243 

0.13243 

0.1324257 

0.13243 

'0*8 75 

0*069854 

0.069854 

0.0698503 

0.069853 

Smaller 

solution 

corresponds to the 

exact solution 

with 


B = 3—2 2, 
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TABLE 6^.5 

Numerical Results for Problem 6*3 

(*■ = 0 , 8 ) 


X 


h 


' " " 

1 

ro • 

i 

5o 

1 

40 

1 

60 

0 ,0 

0.2428968 

0.24012132 

0.2394465 

0.2392385 

0,1 

0.2,3842 25 

0.2 37062 3 

0.2367227 

0.2366591 

0,2 

0*2303030 

0.2293343 

0*2290929 

0*2290477 

0 ,3 

0 .217388 3 

0*2166573 

0.2164754 

0*2164412 

0.4 

0*1996513 

0,1990947 

0*1989563 

0 *1989 3Q2 

0 ,5 

0*1771833 

0,1767655 

0*1766617 

0.1766421 

0*6 

0*1501 282 

0*1498253 

0,1497502 

0*149 7360 

0 *7 

0.1186656 

0.1184597 

0.1184087 

0.1183990 

0 *8 

0 *08 3002 3 

0.0828782 

0*0828474 

0 *0828416 

0 .9 

0*0433660 

0*0433102 

0*0432964 

0*0432938 



Numerical Results for Problem 6 *4 


( (3 = <p = r = l) 


X 

1 

10 

1 

20 

1 

40 

1 

So 

0*0 

0*8 359614 

0*8 365633 

0*8 367104 

0 *8367368 

0*1 

0*8 38-2945 

0 *8 38 346 3 

0-8 383594 

0*8 38 3613 

0 *2 

0*8431136 

0 *8431656 

0*8 4 31787 

0-8431808 

0*3 

0*8511594 

0 *8512115 

0*8512247 

0-6512270 

0 *4 

0*3624522 

0 *8625039 

0 *8625172 

0*8625196 

0 .5 

0 *8 770181 

0 *8770684 

0-3 770817 

0*3770842 

0*6 

0*3948878 

0 *8949 352 

0*8949473 

0*3949505 

0 *7 

0*9 1609 37 

0.9161358 

0*9161470 

0.9161497 

0 *3 

0 .940666 3 

0.9406997 

0*9407037 

0*9407107 

0 .9 

0*9636 306 

0 *9636505 

0*9686559 

0-9636570 
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Table 6 .7 



Numerical 

Results for 

Problem 6 .4 



o = 

0 *05 , r = 20, cp = 6 ) 


X 

1 

10 

1 

To 

h 

1 

40 

1 

60 

0.0 

0 . 0015533 

0.0017349 

0.0017841 

0 .0013013 

0.1 

0 .0013073 

0.0020163 

0.0020326 

0 .00209 4 3 

0.2 

0 *002 7098 

0.0030297 

0.0031299 

0.0031475 

0.3 

0*0048 2 77 

0.0053973 

0 .0055753 

0.0056065 

0.4 

0.009 5 302 

0.0106401 

0.0109880 

0.0110492 

0.5 

0.0200037 

0.0222778 

0 * 02299 32 

0*0231190 

0 .6 

0*0436007 

0.048 3495 

0.0498419 

0*0501037 

0*7 

0.0971114 

0.1068 313 

0.1098473 

0.1103717 

0 *8 

0.2175953 

0*235.6493 

0 .2410369 

0.2419535 

0 .9 

0.4795661 

0.5039730 

0.5106301 

0.5117135 



CHAPTER VII 


IMBEDDING METHODS TO 
INFINITE INTERVAL PROBLEMS 

7*1 Introauction, ; The object of this chapter is to describe the 
a. p p iro a cn of asymptotic bounaary condition to solve boundary 
value problems on semi- in finite intervals- As discussed in 
Chapter 1/ early attempts to treat these problems aro to 
replace them by finite interval problems/ before confuting 
the solutions- Also one can reduce the problem to finite 
interval by coordinate transformation* 

Following the lines of Fox [42] a second order finite 
difference method has been discussed by Robertson [so] for a 
boundary value problem on infinite interval/ by truncating 
it to a finite interval problem- Truncation procedure has 
been applied by many authors* Similarly using the coordinate 
transformation technique/ Grosch and Oraszag [46] presented 
a Chebyshev polynomial approach to solve some physical problems 
on infinite intervals* The author reduces the infinite 
interval problems to finite interval on [— l/l] using the 
transformations and employed the Chebyshev polynomial approach 
to solve the problems* Many authors also followed the 
transformation procedure to solve these problems by ignoring 
the effect of singularity. Recently Lentini and Keller [6 2] 
have analyzed these problems by determining appropriate asymptotic 
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boundary conditions at a finite point. Marbowich [64] has 
also described an adhoc method to solve boundary value problems 
which are posed on infinite intervals by reducing the infinite 
interval to a finite but large and to impose additional boundary 
condition at the far end. 

In this chapter y we have presented an asymptotic 
boundary condition for the numerical solution of the point 
boundary value problem posed on infinite interval. By reducing 
the infinite interval to a finite interval which is large and 
imposing appropriate asymptotic boundary condition at the 
finite point, the resulting boundary value problem is treated 
by using invariant imbedding methods* The tridiagonal system 
resulting from imbedding is effeciently solved. The stability 
of the method is analysed. The transformation technique has 
been briefly included as an alternative method. The numerical 
solutions for the model problems have also been presented. 

7 *2 Asymptotic Boundary Condition : 

We consider boundary value problems on infinite intervals 
in the following way s 

Ly(x) = y"(x)+f(x)y'(x)+g(x)y(x)=r(x) y a < x < ~ 

sub j ect to boundary conditions 

y(a) = a 

y( °°) = 0 

lim y ( x) = 0 
x-* 00 


or 


(7*2) 

(7.3) 
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where the functions f(x), a(v) i u \ 

t gvxj ana r(x) are continuous and 

g(x) < O » We rewrite the Ecrn* f 7 ^ 1 ) , 0 = , , 

qn* w*i ) as a first order system 

in the form : 

Let v( x) - u(x), Y / ( x) = u f (x) = v(x), we have 

u^x) = v(x) ( 7 . 4 ) 

v'(x) = -f(x) v(x) - g(x)u(x)+r(x) ( 7 . 5 ) 

and correspondingly (7*2) and (7»?) become 


u(a) = a ( 7 . 6 ) 

lim u(x) = Uqq = £ {l ml) 

X->°o 


"t 

Letting U = (u,v) , (t denotes transpose), we can write the 
first order system (7*4) - (7*5) in the matrix vector form 


where 

A(x) 


U' = A(x)U + B(x) 
= F(x,u) 


(7.8) 

(7-9) 


r 0 

i 

— f ( x) 




B(x) 




A general theory for linear and nonlinear systems of the form 
(7*8) on semi— inf inite interval has been discussed by Lentini 
and Keller [6 2]* We assume that (i) lim A(x) = A, a constant 

x -*co 

matrix, (ii) lim = 0 , (iii) A(x) is piecewise continuously 

St-*' 00 

differentiable on (a,oo), and u . is required to be .the root of 
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lim F(x,u) o 0. w _ , , 

x -* od we also assume that A is assumed in 

the canonical form such that A = BJlf 1 / o (a zero matrix) and 
J has the block diagonal form J = diag (J + , j°, j~) wher e J + 
contains eigenvalues of A wit h positive real part, J° the 
eigenvalues of A with zero real part and J the eigenvalues 
of A with a negative real part. The main idea is to find all 
bounded solutions and to eliminate the contribution from the 
unbounded solution of the equation (7«8)» The behaviour at 
infinity of the solution of the system (7«8) is essentially 
given by the eigenvalues of the matrix 


A inf = 1±m A (x) 

t '“’ X -* 00 


0 

K 


1 » 

f 

t 



where K = Lt - f ( x) and L = 11m - g(x) • Suppose the matrix A o 
x-»°° x-*» 

has the eigenvalues A^ and then depending upon positive 

real parts of the eigenvalues Real A^, Real A 2 > 0 we 
the linearly independent solution which decay exponentially 


at infinity and the linearly independent solutions which are 
bounded as x -* Since we need only one condition at the 

far end, we expect only one eigenvalue with positive real part 
and say the eigenvalue is A^* Now corresponding to this 
eigenvalue we introduce the projection matrix P M of the form 
P M = [l/0]* ( If the eigenvalue A 2 is with positive real part, 
we choose P M = [ 0, l] ) * Correspondingly the eigenvector to Aj_ 
and A 2 be denoted by Ej. and E 2 respectively such that 
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E 1 = [ E ll' E 2lJ ' and E 2 Then the matrix of 

eigenvectors of A is. 




J 21 


E. 


22 i 


This essentially allows us to find E X , (the inverse of E) , for 
which we know that E^ 1 A ±nf E = diag (A^), we write the 
asymptotic boundary condition as 


lim P M E~ X F( X /u) = 0 (7-10) 

X-*oo 

Eqn-» (7*10) yields the condition at x = N, where N is chosen by 
taking different values of x for which the computed solution 
approximates the actual solution * 

7*3 Invariant Imbedding t 

In order to solve the finite interval problem, we 
extend the discrete invariant imbedding technique which is 
described in Chapter 2- For the sake of brevity, we assume that 
the additional boundary condition for the Eqn* (7*1) with (7*2) 
is of the form 

% y(xj + 3 0 Y'(xJ = k (7.1l) 

where a Q , /3 Q and k are constants, x^ = N, N large tut finite, 
with a o p o > 0 and |0 0 l + t3 0 t A ^ is guarantees the unique 

solution of the two point boundary value problems given by 
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Eqns * (7-1), (7,2) and (7,11) exists- (cf. Keller [55]). 
Discretizing Eqn, (7*l) and using central difference formulae, 

we get after some manipulation a finite difference system of the 
form 

A i y i+i “ B i y i + c i Y i ^ 1 = D ±/ i = 1,?. (7.12) 

where 

A ± = 1 + h f ± /2 
B ± = 2 - h 2 g ± 
c ± = 1 - h f./2 
D i = h2 r± - 

The corresponding boundary conditions (7-2) and (7.11) can be 
written in the discretized form as. 


or 


y Q = a 


% y n + '®o (-^-^-) - X 


2h a o y n + Po y n-fl ~ ^ y n-l = 2hk 


(7.13) 


(7-14) 


where the value of y n+1 at the fictious point x = (n+l)h needs 
to be eliminated between the Eqns. (7*12), (7*13) and (7»14) * 
We seek a difference relation of the form 


y±+i = y± w i + s i- 1 = °' 1 ' 2 0-1 


(7.15) 
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where VL and S 


i co rrQ sponds to W(x ± ) and S(x ± ) and 


determined. Using Eqns , (7*15) and (7-12) 


we have 


y i x , y> + _^i s i ~ D i^ 

(B i ~A i w i ) ' 1 ~ 1 { EL - 


Also from Eqn.(7*15) we get 


Yi = W, 


,~i y^i + s ± 


We have from Eqns- (7.16) and (7-17) we get 


W 


C i 


i~l 


( B ± - A ± W ± ) 


5 . m { h s i:^. 

>1 ~ 1 ( BjL - A ±Wl ) 


To solve these rocurronco relation for and ( 

we need to know the values of W„ . and S , From 

n~l n-1 

we have 


2h v 2na o 

Y n+i = y n - t + r a * ~ To y » 


Writing Eqn. (7 ->12) at x = n ( and eliminating Y n+ ^ 
Eqn* (7*20), we get 


< A n +C n ) 


2hA n k 

Dn) _ 


2A ha 

(-J-- * V 

H o 


V-l + V ' '2A n ha; 
: n + "K 


n "~° • — " ( B„ + -£'*> 


are to be 

(7-16) 

(7.17) 

(7.13) 

(7-19 ) 

i = n~2, » .o/O ) j 
Eqn. (7.14) 

(7.20) 
by using 


^n * 


(7.21) 
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Using Eqn. (7.17) at x - n and comparing with Sjn. (7.21) « 
have 


W. 


( V«n> 


n-1 


2A ha " 

(B n + -A.- £.) 


(7.22) 


2hA k 

'■ sf - ' V 


n-1 


2Aha 

(Bn + 


(7.23) 


Starting with these values of and S n _^, the values of 

and (i = n— 2,n— 3, *•» *,0) are obtained by backward substitution 
using Eqns. (7*18) - (7*19). The solution s (i = 1,2, ..*,n) 
can then be obtained from Eqn. (7*15) by using the known values 
of W^' s and S^/s (i - 0,1,2, *.*, n-1 ) and the known initial 
value y c » 

"S ' 

7*4 Stability : 

We will now show that the method always provides a 
solution and is computationally stable. We examine the 
recurrence relation given by Eqn* (7*18) <* Assuming that a 
small error E ± has been introduced in the calculation of W ± , we 
have 


W i = W ± + E ± 


(7.24) 


and we are actually solving 


w 


c. 

1 


i-1 


(B ± - A ± W ± ) 


(7.25) 
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From Eqns « (7,25) and (7.24) we get 


E, - ~ W. 


W, H E. 


i-1 - i~l cT 

1 x 


(7.26) 


under the assumption that initially the error is small. To 
show I 1 < l E il/ we first show that 0 < < 1 , for 

i = n-1, .,.,0. Assume that A ± > o and C ± > o for all 1 and 
from the definition of A ± , and C ± , we also have B ± > 0 and 
B i > ^i* 


n-1 “ 2A^ha ' 

( -V* + ^ 

we see that under the abore rationed conditions 0 < w „ < 1- 

n-1 

A1 so 


( B 1 - A W_ 1 ) 

n-1 n— 1 n— 1 


< lB n-l " Vl> 


< 1 

and thus 0 < W 2 < 1. Similarly, it can recursively shown 
that 0 < Wj_ < 1, i = n— 3, » » »,0 » Thus the error Egn. (7-26) 
gives | E^ < l E^i , provided A ^ and therefore the 

recurrence relation (7 .22) is computationally stable* Following 
a similar arguments it can be shown that Egn* (7*23) is also 
stable# 
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7 «5 Coordinate Transformation ? 

In order to demonstrate the transformation procedure 
over infinite interval, we have considered different transformation 
for the assumed Egn. (7,1) ~ ( 7 . 3 ), A standard procedure is 
let X = a/x, if a / 0 or X = if a = 0. By using algebraic 

transformation one may reduce Eqn. (7*1) as follows : 


X__ X 


(7.27) 


dX 

dx 


( x+1 ) 


= ~X 


& = „ x 2 dy 

dx dX 


,2 


dx 




dX 


(7.28) 

(7«29) 


By using Eqns. (7*28) - (7»29), Eqn. (7.1) becomes 

+ S i| - Y(X) - s'f- 5 - 


.} + y(x) = 

X X' 


( 7 . 30 ) 


X 


( x) 

where f(x) = £(x), g(x) = g(x) and r(x) = r(x) with x 

The corresponding boundary conditions Eqns. (7*2) and (7»3) becomes 


(7.31) 

(7.32) 


y(x) = a 
y(o) = 0 

where X is defined in Eqn. (7.27) * 

Similarly using an exponential transformation one may 
reduce the infinite interval problems to finite interval problems 
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Using X _ 1-e , Eqn- (7>l) reduces to 



(7-33) 

(7 > 

(7-35) 


Eqns. (7.33) - (7. 35) and Egns . (7*30) - (7-32) constitutes a 
singular two point boundary value problem either with a regular 
singular point or singularity of the second kind. If one is led 
to the regular singular problem, then the methods discussed in 
Chapters 2—3 can be applied directly, by using the expansion 
procedure about the singular point* If the transformed problem 
has singularity of the second kind, one may apply a different 
procedure (cf* de Hoog and Weiss [32]). However, in general, 
these types of transformations produce a singularity of the 


second kind for infinite interval problems.-. Depending upon the 


physical situation one may choose a proper algebraic or 
exponential transformation * An obvious advantage of the change 
of variable approach is that no experimental choice of a large 
number to represent the point at infinity is necessary* This 
particular procedure has been applied to the problem 7*1 and the 
results are tabulated in Table 7-4« 
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7 *6 Numerical examples and Discussion : 

In this section some numerical examples of the model 
problems are solved to show the efficiency of the method* 

Problem 7*1 : ^o illustrate our method/ we solve 

Ly = ~y" - 2y / + 2y = e~ 2X / o < x < ~ 
subject to boundary conditions 

y(0) = l.o 
y(<») =0*0 

This problem has earlier been solved by Robertson [so] and its 

exact solution is given by y(x) = \ e'^ + ^ x + e~~ 2x * The 

asymptotic boundary condition for this problem by our method is 
given by 2y(x co ) - ^3 y'Cx^,) =0* 

Problem 7*2 : As a last experiment, we solve 

~y* + ( 1 + ~) y = , l<x<~ 

with boundary conditions 

y(l) = 0 

y(co) — 0 • 

Pope [42] and Robertson [80] have earlier solved this model 
problem. The asymptotic boundary condition for this problem 
by our method is given by 0*5 yCx^) +0*5 - 0* 
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The numerical results for problem 7.1 for different 
values of N are presented in Table 7 . u ^ values of N 

for computation are 8 , 9 and 10. We observe from this table 
that the computed solutions compare well with the exact solution. 
Also the solutions are seen to compare well upto five to six 
decimal places of accuracy with Robertson's [80] solution at 
Xoo = 9 • The errors present in the numerical solution for 
' these values of N are given in Table 7-2* It can be observed 
from results that = N = 9 can be taken to represent the point 
at infinity for this problem* 

The numerical results for Problem 7*2 are shown in 
Table 7*3 considering the variations in the computed solutions 
at different x for different step sizes, we observe that the 
value of 3 ^ = N gets smaller for smaller values of h* It is 
evident from the table that the solutions decay relatively 
slowly which agrees with the observation made by Pox [42]* The 
computed solutions also compare well with that obtained by 
Robertson [80] at = N = 9* Prom the computational results, it 
thus can be stated that N = 9 represents the point at infinity for 
this problem* 

In Table 7 *4 we present the numerical solutions for the 
Problem 7*1 using the coordinate transformation* We applied the 
exponential transformation to get the finite interval problem 
and used the expansion procedure near to singularity which is 

« The value of X = 6 as used in 


discussed in Chapter 3 
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Chapter 3 f has been taken near' -kwo »-,• 

n ar the singularity (viz. 6 = 0 . 9 ) 

to observe the variation -in 4-v,^ . . 

3tl0n “ the solutions by transformation. 

The computed solutions compare well with the exact solution up 

to five places of decimal. Though the latter method reduces th< 

infinite interval problem to a problem on a finite interval 

thereby evading the computation on a larger interval (as in the 

former case) but the main disadvantage is that it leads, in 

general, to a singular problem with an essential singularity. 

In view of this, the technique of coordinate transformations 

is not practical in many cases. 
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TABLE 7 .2 


Errors 


in the numerical solution for Problem 7-1 
for different values of N 



■ **- ■ - — 

• - - 

- - - 

Xqo = N X 

h 

1/64 

1/128 

2.0 

0 #165 E— 05 

0.171 

E-06 

4.0 

8 

0.29 3 E-07 

0.325 

E-08 

6 .0 

0.323 E-08 

0.279 

E-08 

o 

* 

CO 

0.119 E-07 

0.119 

E-08 

2.0 

0.165 E-05 

0.171 

E-06 

4 *0 

0.287 E-07 

0.264 

E-09 

9 6 .0 

0.661 E-09 

0*219 

E-09 

a ,o 

0.78 2 E-09 

0.774 

E-09 

9 .0 

0.161 E-09 

0.161 

E-09 
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CHAPTER VIII 


CUbi pr^M 1 im? Solutions of Boundary Value 
Problems over Infinite Interval 

8 .1 Introduction : The development of improved methods for 
solving two point boundary value problems of practical 
significance is of great important in Science and Engineering. 
In many cases the domain of the governing equations of these 
problems is infinite or semi- infinite so that the special 
treatment is required to these so called infinite interval 


problems# These problems occur very frequently and are of 
groat importance in areas such as fluid dynamics,, aerodynamics 
etcetera# Often in most cases the analytical solutions for 
these problems are not readily attainable and thus the problem 
is brought to the problem of finding efficient computational 
algorithms for obtaining the numerical solution. 


Ever since the pioneering work by Ahlberg et al* [ 2 ] 
there has been a great deal of development in the theory of 
spline functions and their applications to several practical 
problems* Of these, the cubic splines have attained a prime 
place and have attracted the attention of many to solve, in 
particular, boundary value problems# This is possibly because 
those cubic splines are efficient and simple to use and possess 
important properties that are required of a good approximation# 
To cite a few, Bickley [ H] has considered the use of cubic 
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spline for solving linear two 


point boundary value problems. 


which leads to the solution of a 


set of linear equations whose 


coefficients matrix is of upper Hessenberg form. The cubic 
spline method suggested by Bickley has also been examined by 
Pyfe [ 45 ] in combination with deferred correction to solve 
two point boundary value problems. The cubic spline approximation 
has been applied by Albaslny and Hoskins [ 4 ] to an integral 
equation reformulation of the original differential equation, 
which have been shown to have smaller truncation error than 
would be obtained by direct use of the cubic spline on the 
differential, equation itself* In another paper by the same 
authors [3], the spline solutions have been obtained by solving 
a set of equation with a tridiagonal matrix of coeff icients - 
Daniel [28] has proposed the use of acceleration with collocation 
with various types of spline approximation appropriate for two 
point boundary value problems* There are several other papers 
on this topic, but little seems to have been done in using 
cubic splines to solve boundary value problems over infinite 
intervals • 


In "this chapter we present a cubic spline procedure for 
numerically solving boundary value problems over infinite 
intervals. Unlike the familiar three point finite difference 
scheme second order accuracy is maintained even with a non- 
uniform mesh and relatively large changes in the grid spacing. 
Also spline approximations described here leads to tridiagonal 
systems *» The asymptotic boundary condition for the infinite 
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interval problems is derived mu 

derived. The cubic spline formulation and 

the procedure is discussed with derive™ w. 

derivative boundary conditions. 

The tridiagonal system has been solved efficiently and the 
stability of the method is discussed. Numerical examples for 


the model problems are pivon 

8 .2 Asymptotic Boundary Condition : 

We consider the linear two point boundary value problems 
of the form 


Ly(x) = y*(x)+P(x)y / (x)+Q(x)y(x) = R(x), a < x < ~ (8.l) 

subject to boundary conditions 


y(a) = a 
y(°°) = j3 

or lim y(x) = 8 
x-*°° 


( 8 . 2 ) 


(S.3) 


where the function P(x), Q(x) and R(x) are continuous and 
Q(x) <0* In order to find the appropriate boundary conditions 
for the Eqn. (8*l) we rewrite (8*l) as a first order system 
in the form 

u'(x) = v(x) (8.4) 

v'(x) = -P(x)v(x)-Q(x)u(x)+R(x) (8*5) 

where y(x) <= u(x) J y' (x) = u # (x) = v(x)« Correspondingly ( 8 . 2 ) 
(8*3) became 
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u'(a) = a ( 8 , 6 ) 

Lt u(x) = 0 ( 8 , 7 ) 

X-*co 

Letting U = (u,v) , (t denotes transpose), we can write the first 
order system (8*4) — (8 #5) in the matrix vector form 


U' = A(x)U + f(x) 


(3.8) 


where 


= F(x,U) 


A(x) 


j i f(x) = l | 

j^-Q(x) ~P(x)! , ] R(x)j 


As discussed in the earlier chapter, the method requires that 
the matrix A(x) should be a constant matrix as x •* «>« Also the 
main idea is to eliminate contributions from the unbounded solution 
of the equation, which are essentially given by the eigenvalues of 
the matrix A(x) * Let denote the matrix when x -* 00 so that 


A 


inf 


Lt A(x) 

X-+0O 


where K = Lt - Q(x) and L 

X-*oo 


\ 0 1 | 

| K L j 

L. — ► 

= Lt - P(x) • 

x -*co 


Suppose the matrix 


A inf has the eigenvalues A^ and .A 2 , then depending upon Real X^ 
Real A 0 > 0, one can find linearly independent solutions at 
infinity and the corresponding projection matrix. Let E be the 
matrix of eigenvectors of A ±nf/ then by calculating E" 1 for which 
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E A inf E “ dia 9 we write tiie asymptotic condition as 

lim P M E" 1 F(x,U) =0 (8*9) 

X-*oo 

: Eqn. (8*9) yields the condition at x = N where H is chosen by 
talcing different values of x for which the computed solution 
approximates the actual solution. 

8 *3 Cubic Spline Formulation : 

We now describe a cubic spline procedure to solve the 
finite interval problem obtained by asymptotic boundary condition* 
For the sake of brevity,, we assume that the Eqn* (8*9) is of the 
form 


a o + £ 0 = K 


( 8 * 10 ) 


for x^, = N, N large but finite, where a Q , $ Q and K are known 
constants such that <2 O 0 O ° and 1<X 0 I + 18 0 I 0* Ihis guarantees 

the unique solution of the two point boundary value problem 
given by Eqns* (8*l), (8.2) and <8*10) (cf. Keller [ 55 ]) * 


We consider a mesh with grid points a = < x^ < x 2 *»»<x n =N 

with h = ^ > 0* The function y(x) at the node points x^ 

is denoted by y(x^) = yj_/ Let S Q (y/x) = S^(x) denote the cubic 
spline function, which is continuous together with its first and 
second derivatives on £a,N] and satisfies S^(x.^) = Yi* Then 
in general : 


SJ( 


X 


) = M 


i-1 


( X£~x) 

I ‘ "h ' 


“i 


( X “ X i-l } 

K * 


( 8 . 11 ) 


138 


where M i = S^(x i )* 

Integrating Eqn* (8.1l) twice results in the spline 
interpolation formula on (x^/Xj^) t 


( x^-tx) J ( x^x ^ - ) 

S p (x) = M j__i 5‘h + M i 6h 


2 ' M 4 _ 4 h^_ (x^"x) 


• + tyi-j 


I"! 

6 ; h 


M.h“ (x— x. 1 ) 

♦ ( y ± - -f” 5 -~ i h - 1 ” 


(8-12) 


where the integration constants have been evaluated by the 

requirement of continuity of the spline function and its first 

derivative at the node points* Ahlberg et al» [ 2 ] have shown 

that if the function y(x) e C 4 [a,N], then the spline function 

S ( x) approximates v(x) at all points in [a.*N] to fourth order 
3? 

in h# The unknown derivatives are related by enforcing the 
continuity condition on S p ( x) » 

Differentiating (8*12)., we get 


(x.j-x} 2 h (x-x. .,) 2 hi ^ y i” y i-l^ 

Sp(x) - M ±-1 [- — 2h + F 3 + M i^ " 2h " + ' h ’ 

(8*13) 

From Eqn* (8*13) we have the one sided limits of the derivatives 


as 


5 '(*f) = - 1 % " F M l + 1' 1 " °' 1 ' 2 ' 

Jr 


(8*14) 


and 


S'(^) : 


( yi _ Y irli. + h . + ^M. iy i = 1/2* 

3 i 3 i*-* X 


.**n (8-15) 
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By virtue of (8*ll) and (8*12) the functions S*(x) and 
Sp(x) are continuous on [a,Nj, the continuity of S|^(x) at x^ 
yields by means of Eqns. (8.14) - (8.15) the condition as 


l M i-i 


20 ! M . + *2 M . 

+ 3 x + BT u i + i 


y i+r y i 


h ^ 


i = 1/ 2* » *#,n~l 
(8.16) 


Employing spline approximation (8.13) to the Eqn. (8#l) and using 
(8 *14) - (8.15)/ we have for i = 0,1/2/ ..^n-1/ 


»i + Pi [-If- - 1 - £ M i+1 ] + Q i y i = R ± 

or 

(i- | PjM 1+1 + (y 1+1 -yi)«iyi = R ± (3 - 17) 

and similarly for i = 1,2, «.*,n, 

P 

(i + |p 1 )m 1 + § + if ( yr y i-i )+Q i y i = R i (8 * 18) 

Also from Eqns* (3«17) and (8*18)/ we get a relationship 


2 Mi + £ P i M ± _ 1 


5 P l M i + l + 


h 


'yi+i-yi-i i * * * 5 = 2R i- 2Q i y i 


i ss I, 2/ * * * , n~l • ( 8 *19 ) 

To find the cubic polynomial S^(x) at different points we need 

the unknown coefficients and M ±+1 in the Eqn. (8*12). 

We eliminate M^s ( i = 0,1, . ..,n) from Eqns. (8.17) and (8.18). 

From Eqn. (8.17), we replace i by i-1 and get 
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a -§ p i-i ,M i-i-F p ±_i M i + T 4 = Vi 


( 8 . 20 ) 


From Eqns. ( 8 » 20 ) and (8.18) by eliminating we have,, 
M i~l = FT ^ P ±-i + l P i R i-l P i-l R i ) 


, P P.P. 

*■* y^^ ( 2 , 3 P i^i~l ”* * H' *" 2 1 ') 


P P P 

— y^( ^ ^ “** ^ ^ -1 /^‘ == ^ / # * #/n (8#2l) 


Where a ± = 1 ~ § P.^ + f P i ~ P i P i-l 


( 8 . 22 ) 


Again from Eqn. (8-17), we replace i by (i+l) and get. 


l i+l = FT £ (P i+l-| P i R i+l"F P i+l R i ) 


b. r.. + ^ - ™o 


y i+l (Q i+l ” 3 P i R i +1 + ‘ h 


P^P. 


+ y^( * + ^ P i+l®i^ * i= 0 /l/ 2 , 


(8 .23) 


with b ± = 1 - 3 Pj_ + 3 P ±+1 - p i p i+i 


(8.24) 


By substracting Eqn. (8.16) from Eqn. (8.19) we get an expression 
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y i~l (l ~ 1 P i ) - y i (2 - Q i } + Y i+l (l + f p i ) 


2h ’ R, + 4- (1 - 3 P i )M i-i + T" (1 + 


3 i T 6 


3 i i+1 


(8*25) 


i = lr 2/ * * * # n—l * 

Using Eqns • (8*2l) and (8*23) in (8.25), we have after some 
algebric manipulation a tridoagonal system given by 


D i y i-l C i y i + E i Y i+1 = F i / ^ 2 ' * * (8*26) 

where D ±/ C . T and are given by the equations 


= b^ [l — ^ P i—1 "** ~"S" ^i— i ^ -i 


(8 *27) 


i = C a i^l + \ p i+i^ +b i (1 - 2 P i~l^ 


and 


2h 2 n h£ p P + 7h ( P P . ) ] (8.28) 

~ ' 3 ' Q i U 12 ^i-l i+1 + 24 i+1 i-l 


E 


i = a i C 1 + $ P i+1 + IT Q i+1^ 


h z 


(8.29) 


E ± = | [ h 2 b^R ;L ^+4h 2 R^( 1 - ^ P i~l P i+l + 24 (P i+l” P l-*l^ 


+ h a i R i+1 ] 


(8.30V t 


*‘V 


8 .4 Computational Method : 

To solve the tridiagonal system given by (8.26) we seek 
a difference relation of the form 
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— W^y^ + i — 1 , 2, • «-./n— 1 (8 ■•31) 

where and correspond to W(x^) .and T(x^) and are to be 
determined* By using (8*3l) in (8.26) # we have 




( C i- E i w i> 




- E ±Wl ) 


But from (8.31)/ we have 


(8.32) 


Y± = W i~l Y i-1 + T i-1' 

From Eqns. (8-32) and (8.33), we get 

D. 

1 

w. = — • — 

1-1 ( c i - E i w i.) 

and 

< E i?i - p i? 

T ss - - 1 

i " 1 (C ± - EjW ± ) 

To solve the recurrence relation for W i and 
we need to know the values of W n-1 and T n-1 ' 

From Egn. (8.10) for x = n, we have 


(8.33) 


(8.34) 


(8.35) 

( i = n— 2, . * . , 0 5 


[ n 



(8 .3d J 


Eqn. (8.36) can be approximated at x = n by using the expressions 
(8.15), (8.21) and (8.23), we have 

a i Yn ~ ^1 y n-l = r i 


(8.37) 


143 


a two term relationship in y n and y n , where 


“1 - «o + 3a n Po<-°n ~ TT ~ ' Sr + t P n-1 °n + ’ V* * v!> 
LX n 

(8.38) 


B - h ’ b (- fs - Pn “' 1 + °n-l P n P n~l 3 a n% 
" 3a_ *V“ ~ 2h + ~"2 * + “4 “ + ; r 0 


n 


h 


(8.39) 


and 


T 1 - K + Jt ^-Xn Vl + X P n-l R n ) 
n . 5 - 


(8.40) 


From (8.33) for i as n we have 


y_=W - Y - +T 
■*n n-1 ^n-l n-1 


(8.41) 


Comparing (8.37) and (8.41),, we obtain and as 


W 


n— 1 = ^ a l 
‘n-1 = T l/ Ct l 


(8.42) 

(8.43) 


Then W^s and T/s (i = n-2,n-3, •• *,0) are obtained recursively 
in the backward process by using Eqns. (8.34) — (8.35). Using 
the values of W/ s and T^s and knowing the value of y D , the 
solution y^/s (i = 0,1,. ..,n) can be obtained by forward process 
by using (8*3l). 

8.5 Stability : We now examine the recurrence relation given 
by Eqns« (8 . 34) - (8.35) for stability. Assume that a small 
error e ± has been introduced in the calculation of W ± then we 


have 
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W i = w i + e i 
and we are actually solving 

«s> ^ . 

W. . = - ■ - 1 - V- 

' c i - E i «±> 

From Eqns • (8 *34) and (8-45), we have 


6 i-l “ W i-1 e i W i-1 7 

under the assumption that initially the error is small. 


(8-44) 


(8-45) 


(8- -46) 


Let us assume that E^ > 0 and > 0 for 1 < i < n~l, then 
from the definition of Djy and since Q(x) < 0, it can be 

easily verified that > Dj+E^ for 1 < i < n-~l • From Eqn- (8-42) 
we have W n ^ = 8^/a^., and < 1 , if M > 0 and 8^ > -M/2, 

where M = a Q + 0 Q (~Q n - + § P n _ 1 Q n ) - Under this 

Ja n 

condition and making use of the assumptions on D^ f and E^, 
it follows very easily from (8*34) that 


I Wf| <1/ for i *s n— 2,n~3, » • *,0 

and thus (W^i <1/ i = 0/1, 2, • » -,n-- 1 * From Eqn- (8*46), it 
then follows that 

2 1 E i- 

^i-l* = !W i-l ! ' * ,e i' 

11 1 A i Dj^| 

< provided [E^t < 1°^! 


(8-47) 
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making the recurrence relation (8 #34) stable# Similar arguments 
will show that the recurrence relation (8 #35) is also stable# 

8 #6 Numerical examples and Discussion s 

In this section, the numerical results of the two 
model problems solved are given to illustrate the method and 
its effectiveness# 

Problem 8*1 : We solve 

Ly( x) = -y" - + 2y = e~^ x , 0 < x < 

with boundary conditions 

y(0) =1 

Lt yCx) 0* 
x-«> 

This problem has earlier been solved by Robertson [80 ] and its 
exact solution is given byy(x) =2 e + ^ e The 

asymptotic condition obtained by our method for this example is 
given by 

-r ytaO + (-1 + ^3)^ (aq*) = o# 

^ 2V3 

Problem 8 #2 : As a second example, we solve 

Ly( x) = — + (1 + ~)y - -^2 j 1 < x < °° 

X 

subject to boundary conditions 

y(l) = 0 

lim y( x) =0# 
x-»°° 
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This problem has earlier been considered by Pox [42] and later 
by Robertson [80 ]• The asymptotic boundary condition is given 
by \ y(xj * \ y'(x«,) = 0, 

Table 8*1 represents numerical solutions for Problem 8*1 
at some selected points by taking different mesh size h and 
different values of = N* The values of N taken for 
computation are N = 10, 11 and 12* The exact solution of the 
problem is also presented and it is observed that the computed 
solutions compare favourably well with the exact solution* It 
is known for this example that the solution decays exponentially 
and this fact is reflected in the computed solutions -for 
different values of N* It can be seen that the computed solutions 
for N = 12 give atleast six place occuracy and thus N = 12 can 
be taken to represent the point at infinity for this problem* 

The numerical solutions for Problem 8*2 for different 
values of h are given in Table 8*2. It is observed from this 
table that the solutions decay relatively slowly which agrees 
with the observation made by Pox [42 ] * The computed solutions 
also compare very well with that obtained by Robertson [80 ] and 
at x^ = N =9, the solution compare upto five places of 
decimal* The advantage of using cubic spline is that it not 
only gives an approximation to the solution y(x) but also an 
approximation to the derivative y'(x) at every point of the 
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LIST OF PROGRAM FOR CHAPTER — XXX 



APPENDIX - III 
************** 


MAI»» PROGRAM 

CDi'-TlMUl)US IMBEDDING FOR SINGULAR 

boundary vai.de problems. 

**************************************************** 

V "m+F! (X)*Y'CX)+F2(X)*YfX)=F3(X) 
YCXo)=ALPHA ; YCbJsBETA . 

*+****+*♦*♦******************************************* 

FI# F2 AND F3 — THE COEFF. FUNCTIONS FROM THE GIVEN EQN 
ALPHA AND BETA — GIVEN ARBlTARY CONSTANTS 
X - CORRESPONDS THE VALUE OF DELTA. 

Ml AND M2 -- ROOTS OF THE JUDICIAL EOUATTOF. 

N1 » HUMBER OF THE TERMS IN THE SERIES EXPANSION. 

Alt., Al(J), 1 = 1,2 Ml AND F10,B1CI) 1 = 1, 2. ..Ml 

C1U,CU.U , 1 = 1,2... .02 AND PJ0,D1<I), 1 = 1, 2. ..Nl 
Af F THE CONSTANTS TH THE SERIES SOLUTION. 

R 1 , R 2 , D R 1 AMD DR2 — CALCULATED BY USING THE 

SERIES SOLUTION. 

DPI AND DP. 2 — THE FIRST AND SECOND DERIVATIVES 

OF THE SERIES EXPANSIONS. 

A OLE = A , BOLE sR AND CONE =C — COEFF.FROM THE BOUNDARY 
CONDITION AT DELTA 


EXTERNAL ALL 

HE A L X , Y C 4 ) , SHUT , RFLERR , ABSERR , K 1 1 , K22 , K3 3 , K44 , 
1 K55 , K6b ,K77 

REAL ZFI NAL , ZPRINT , WORK (39) 

INTEGER I WORK ( 5 ) , IFLAG , NEON 



D f MEHSlUN YK100, 100),U(1 00 ) , V ( 1 GO ) , A 1 ( 1 00 ) , B1 ( 1 00) 
DJMEMSIUW Cl (100) ,I>1 (100) 


RF.AL Ml, ,12 

0PENCUNI1=13,DFVICFs*DSK* , riLE='SUMl .DAT*) 
OF. AD (13,*)N1 ,M1 , M 2 , X 
PEA0(13,*)M0,(AJ(I), T = 1 ,N1 ) 

RKAIH13,*) ,B10, (PKT) , 1 = 1 ,01 ) 
HKAP(13,*)C10,(C1 (I), r = l,M) 

P*: An (l 3,*) DIO, (1)1 (I) , T = 1 ,W1) 

SI 1 — A 1 (F I ) 

DO U I St, HI 
i«2sNj+l-l 
Sn=Slt*XtAl (02) 

sii=sh*x+aio 

PSllsBl(M) 

DO 12 X=?,M 
M 3=01+1-1 

DS1 1=DS1 1 ♦X+Bl CN3) 

DSllsDSl 1 *X+D10 

S22sCl(Nl) 

DO 14 1=2, N1 

M4sNt+l-I 

S?2=S22*X+C1(N4) 

S22=S22*X+C10 



% 


Pane 

r>S22=Dl(iU ) 

DO 15 1=2 , M 
HS~K 1 + 1-1 

DS22-0S22*X+D1 ( N5) 

DS22 = DS22*X + DJ.O 
TYPE*,S11 ,DS 11,522, DS2? 

Hl=X**Mj *SU 

DPl*M!*X**f Ml-1 )*S11+X**M1*DS11 

F2=X**M?*S22 

PP2=M2*X**C02-1 ) *S2?+X**M2*DS22 

TYPE*, HI ,DR1 .R2.PR2 

ADNF=DR2 

BONE=-R? 

CawH*Rl*DK2-DRJ *F2 
TYPE*, A 0 'ME, BONE, CONE 


SAMPLE PROBLEM : 3. 1 

Y(X)=UCX)? Y'tX)U'(X)=V(X). 


Z -CORRESPONDS TO X WHICH IS DELTA TO BE GIVEN, 
Z=C.5 
J = 1 

NUMBER OF EQUATIONS AND INITIAL CONDITIONS. 

NE0N=6 

Y C 1 ) -0 , 0 

YC2)=1 .0 

Y ( 3 3*0 ,0 

Y ( 4 ) = 1 , 0 



Pane 4 


Y(S)=ft.O 

Yf6)=C,rt 


HF%i,EPR= J ,0J>5 
APSEKK=0„n 
ZFTHAL=1 i<j 

X pjjthtso. J 
IKLA«=! 

Zf!JT = Z 

C A Mi RKK4* ( ALL , LEON $Y ,?■ , ZOUT , PELEPR , ABSERR , 

1 IFLAG,WORK, IWORK) 

TYPE* , J , Z , Y ( 1 ) ,Y(2>,¥C3),Y(4),YC5),Y(6) 
on TO c 80 ,26,30,40,50,60,70,80), IFL.AG 
ZOUTrZ+ZPRINT 

on 111 I=),MEON 

Y1CJ,7)=YCX) 

JsJ + 1 

IECZ.r.T.ZFIMADGO TO 10 

SOLVE THE SYSTEM OF SIMULTANUOUS EQUATIONS, 

H11=UL-Y(3) 

VDFLTA= (CM 1 *Y f 4) +Y( 1 ) * Y ( 6 ) ) *AON E-CONE* (Y£ 4)*Y(2)-Y(5)*Y(15)) / 
1 CYC1 )*A0NF-C(Y(4)*Y(?)-Y(5)*YC1))*B0NE)) 

UDFLTA = (CONE-BONE* VDFT.T A) / A ONE 

type* , i/delta * vdelta , vdelta 


U(J )=UOELTA 
VC 1 JaVnELTA 
U (J-l ) »UL 

SOLUTION IN TERMS OF UCD'S. 
DO 551 1=2 , J-2 



Paae 


V( f ) = ( VHKIjT A-Y 1 (T,5)*UDELTA«Yl(l r 6))/n<i,4) 

U( 1 ) = Y 1 (I, n*V(I)+YlCl # 2)*UDFLTA+Yia # 3} 

TYPf* , ll ( J ) , y ( T ) 

CONTINUE 
TYPE 232 

FORMAT ( 5 X , ' — ......... ... • / j 

TYPF 235 

FORMAT ( BX , * INVARIANT IMBEDDING SOLUTION'/) 

TYPE 23A 

FORMAT ( 5X , * — — — ...... .... ... « / > 

DO 555 1 = ) M J-1 

TYPE 233 » l- £ I ) , VI I ) 

format (fix, E20 . 1 A, 5X , E?0. 1 0/ ) 

CONTINUE 
TYPE 237 

FOR MAT £ 5 X , — — ............ » i 

STOP 

TYPE 31 , HELblPR , ABSERF 
GO TO If 
GO TO 1A 
ARSBPRsl . OF-1 3 
TYPE 31 , RFLERR , APSPRP 
GO TO 10 

rflerrsio,o*re:lepp 

TYPE 31 , PEI.ERR , APSFRR 
IFLAGs? 

GO TO 10 
TYPE 71 
IFLAG*2 
GO TO 10 



i-'aae <■» 


TYPE 81 

F(iR'>iAT(2X , ' TOLERANCES RESET ' ,5X,2EJ2.3) 
FORMAT (5X, ' PiUCM OUTPUT * J 
FORMAT C *> X # 'IMPROPER CALL') 

STOP 

EAT) 



SUBROUTINES ARE ALL , PKF45 , RKFS, FFHL:- 


RKF45 - Fehlt-era Fourth-Fifth Order 
Runae-Kutta Method, 

REFER KfJCES : 

E , FEULBERO , J.nw-fiRDFR CLASSICAL PURGE 
KUTTA FORMULA STFPSIZF CONTROL, NASA 
TP-R315. 

I..F.SHAMP1PE, H, A. WATTS. 

SAUDI a LABORATORIES REPORT SAND 7fo-05R5. 


FOR THE TEST EXAMPLE 3.1 THE CORRESPONDING INITIAL 
VALUE PROBLEMS ARE GIVEN . 

YP(1) — Sl'(X) ,YPC?) — S2 ' (X) , YP( 3 ) — S3'(X), 
YPC4) -- Ql'(X), YP(5 ) — 02 ' (X) , YP(6) — Q3'(X). 


SUBROUTINE ALL ( Z , Y , YP) 

REAL Z,Y(4) ,YP(4) 

YP(l)sl,0+(0,5/Z)*Y(l)-(Y(l)*Y(l)) 

YPC2)»-Y(1)*Y(2) 

YP(3)=-Y(1)*Y(3) 

YPC4)=((0.5/Z)-YC1))*Y(4) 

YP(5)=-Y(4)»Y(2) 



Pane 


Y P l 6 ) s -Tc ( 3 ) * Y ( 4 ) 

RETURN 

Eiin 


SUBROUTINE RKF45(F,NF0U ,Y , Z , ZOUT , RELERR , ABSEPR , IFLAG , 

1 WORK, I WORK) 

INTEGER NEON » I EL AG# I WORK (5 ) 

REAL YCNEON) ,Z, ZQUT,RELERR, ABSERR, WORKC1) 

EXTERNAL V 

KIMbNEON-U 

KIbKIM+1 

K?=Kl+NEON 

K3=K2+NEOF 

K4=K3+NEON 

K5-K4+NE0E 

K0-K5+NEON 

CALL RKFS CF # NEON , Y , Z, ZnUT , RELEHR , ABSERR , IFLAG# WORK ( 1) , 
lVfURK(KIM) # WORKCKl),EORK(K2) #WORK(K3) ,WORK(K4) ,WORK(K5) , 

2 WORK (Kf> ) , WORK (K6+3 ) , TV’ORK ( t ) # I WORK ( 2 ) , I WORK ( 3 ) # IWORK (4 ) , 
31WOPK15)) 

RETURN 

EK’D 


SUBROUTINE RKFS ( F , NEON #Y , Z , ZOUT , RELERR # ABSERR , 3 FLAG, YP ,H , 

4F1 # F2 # F3 # F4 ,F5 , SAVRE , tSAVAE ,NFF# KOP # TNTT, JFLAG > KFLAG) 

LOGICAL UFA I LD. OUTPUT . 

INTEGER LEON # TFLAG # NFE # KOP, INIT , , IFLAG , KFLAG 
REA L Y (NEON ),Z. ZOUT , RELF.PR , ABSERR , H , YP ( NEON ) 

REAL F1(NEQN).F2< BE OR )#F3( NEON 3#F4( NEON ) #F5( NEON ) , SAVRE , SAVAE 
REA L A # AE , LT # EE # FEOET , KSTTOL # FT # HM I N # REM IK), RER , S , SCALE , 



Pani 


1 rOL, Tt ILL t tj 2 6 t F PSP 1 f FP,s , Y pK 
EXTERNAL F 

INTEGER K,MAXNFE,MFLAG 
deal AhAXl # AMINl 
data RENIN/3 .r>fi>6/ 

DATA MAXFNF/30OO/ 

IFCNEON.lt. 1 )Gn TO 10 

If C CRLLERP.LT. 0.0) , OP . C ABSERR.LT.O .0 3 ) GO TO 10 
MFLAGsIABSCIFLAG) 

If C (MfTiAG.EO.O) .OR. (HFLAG.GT.8) )G0 TO 10 

IKCMFLAG.NE.l )G0 TO 70 

EPyrl.O 

FtflsKPS/'/.C 

EFSPlsFPS+1 .0 

IFCEPSIM .GT.l) GU TO 5 

D26»2#>.0*EP{S 

GO TO 5P 

IFLAGsft 

RETURN 

IFC CZ.FO.ZOUT) .AND, CKFLAG.NE.3) )G0 TO 10 
IF C FFLAG . NE , ? ) GO TO 25 
IFC CKFLAG.EQ.33 .Op. CINJT.EQ.O))GQ TO 45 
IFCKFLAG.EQ.43GU TO 40 

IFC CKFLAG.E’Q.5) .AND. CARSERF.EQ.Q.O) 3G0 TO 30 
IFC CKFLAG.E0.6) . AND . CRFLERR .LE ,SA VFE) . AND. 

1 CABSERR.LE.5AVAE))G0 TO 30 
GO TO 50 

IFC I FLAG, EQ. 3 3 GO TO 45 
IFCIFLAG.EQ.43GO TO 40 

IF( ( IFLAG.EQ.5) .AND. CABSFRR.GT.0.03 3 GO TO 45 



Pane 


STOP 
NFE = U 

IFCMFliAr. # EP.?)r.O TO 5 0 
IFt.AG=JFLAG 

IF(KFLAO„Ey. 3 )MFLA0sI ABS (IFLAG) 

JFLAGsIFLAG 

KFLAG=d 

SAVREsRFI.KPF 

SAVAEsABSERF 

FER=2 ,0*EPS+RFMIN 

I F ( RELERR , GE • PER ) GO TO 55 

IFLAG=3 

KFLAGs3 

PFTURN 

0Ts(Zni!T'-Z) 

IF ( OFLAG. FQ. 1 )GO TO 60 
1 F ( I M T . EQ • 0 ) GO TO 05 

c;n to be 

iKITst 

KHPssO 

A=Z 

CALL F(AfYfYP) 

HFE=1 

IFCZ.riE.ZOHT)Gn TO 65 

IFLAG=2 

R FT URN 

INTTsl 

Ha(ABSCPT)) 

TOT,fi=0. 

DO 70 Ksl,NEON 



TOL=FELFRR + Ar<SCY (K) 3+ABSERR 
IFCTOL.EU.OJGD TO 70 
TULH=TOL 
YFK=ABSCYPCK) ) 

IFCYPK*H**5.GT.TQL)H=CT0L/YPK)**0.2 
COLT I HUE 

IFC TOLN . LE.0 . 0)H=0 . 0 

H=AMAX1 (H,U?6*AMAX1 C&BS(Z),ABSO)T))) 
JFLAG=ISIGN (2 , JFLAG) 

HsSIGN(H,DT) 

IFCABSCH3.GE .2.0*ABSCDT3)KUP=KOp+l 

IFC KOP.NF . 100 ) GO TO R5 

KOPsO 

TFLAG=7 

RETURN 

IFCABSCDT3 ,GT.U26*ABSCZ)3CO TO 95 
i)D 90 K=1 , HFOH 
YCK3=YCK3+PT*YPCK) 

A=ZOUT 

CALL FC A , Y* YP 3 
MFF=HFE+1 
GO TO 300 
OUTPUT=. FALSE. 

SC A LE=2 . O/RELEFR 

AE=SCALF*ABSEFR 

UFATLD=.FALSF. 

H M I N s U 2 6 * ABS ( Z 3 
DT=ZOUT-Z 

IFCABS(DT3.GE.2.0*ABSCH3)G0 TO 200 
IFCABSCDT3.GT.ABSCH3)GO TO J50 



Pane 


OUTPUTS. TRUE. 

H=|)T 

GO TO 200 
HsG.5*DT 

IF ( NFE , LE . MAXFNE ) GO TO 220 

IFLAG*4 

KFLAG=4 

RETURN 

CALL FEHL(F,NE0N,Y,Z,H,YP,F1,F2,F3,F4,F5,F1) 

NF£»NFE+5 

EEOETsO.O 

DC 250 K=1,NF0N 

ET= ABS ( Y C K ) ) + APS C F l ( K ) ) +AE 

I F ( ET . GT « 0 . 0 ) GO TO 240 

IFLAG=5 

RETURN 

EEsABS((-?090.O*YP(h) + C2197G.0*F3(K)-15048.0*F4(K);>) + 
1 (22528. O*F2(k)-273bO.0*F5(K))) 

EE0F.T=AMAX1 (. EEOET , EE/ET ) 
ERTTOL=ABS(H)*EEOET*SCAI«E/752400.0 
IF ( RSTTOL « LF . 1 . 0 ) GO TO 260 
UFA IIjD*. TRUE. • 

OUTPUT^. FALSE. 

S=C . 1 

IF CESTTOL.LT, 59049. 0)Ss0.9/ESTTDL**0. 2 
H=S*H 

IF(ASSCH) .GT.HMIM)GO TO 200 

IFLAG=6 

KFLAG=6 

RETURN 



Zs7.+ H 

I'O 270 K=1 f NEON 
Y(K)=F1(K) 

A=Z 

CALL F(A,Y,YP) 

NFEsNFE+1 

S-5 .0 

lF<ESTTOL.GT.1.8895e>flE-4DS=0.9/ESTTUL**0.2 

IF(HFAILD)S*AWINi(S,1.0) 

H=SIGN C AHAX1 (S*ABS(H) «HMIN) ,H) 

IF(UUTPIIT)CO TO 300 
IFdFLAC.GT.ODGO TO 100 
IFLAG=~2 
RETURN ■ 

7“Z0HT 
I FL AG=2 
RETURN 
ELD 

SUBROUTINE FEHL(F,NE0N,Y,Z,H,YP,Fi,F2,F3,F4,F5,S) 
INTEGER NEON 

REAL YCNEQN) ,Z,H,YP (NEON), FI (NEON), F2( NEON) , 

IF 3(WE0N) ,F4(NE0N) ,F5CNE0N) ,S(NEQH) 

RIAL CH 

INTEGER K 

CHsH/4.0 

DO 221 K=1,EQN 

FS(K)=Y(K)+CH*YP£K) 

CALF' F(Z+CH,F5,Fl) 

CUs3 .0*H/32 .0 
DO 222 K=1 , NEON 



Pane 


F5CK)=Y(K)+CH*CYP(K) + 3.0*F1 CK)) 

CALL F(Z+3.0*H/8.0,F5,F2) 

CHsH/2 197,0 
DO 223 K=1,NE0F 

F5(K)=Y(K)+CH*(1932.0*YP(fO+(7296.0*F2(K)-7200.0*Ft CK))) 
CALL F(Z+12.0*H/13.0,F5,F3) 

CHsH/4104.0 
DO 224 K-l ,NEON 

F5(K)sY(K)+CH*( (8341.0*YP(K)-845,0*F3(K)) +(29440. 0*F2(K)- 
132832. 0*F1(K)>) 

CALL F(Z+H,F5,F4) 

CHsH/20520,0 
DO 225 jK = l , NEON 

FI (n=Y(K)+CH*f (-6090. 0*YP(K)+( 9295, 0*F3(K) -5643. 0*F4(K})) 
1 + C 4 1 0 40 . n*Fl C K ) -28 352 . C *F2 (K ) ) ) 

CALLF (Z+H/2.L.F1.F5) 

CHsH/761 8050 . n 
DO 230 K = i , titJOL 

R(K)=Y(K)+CH*( ( 902880. 0*YP(K) +( 3855735. 0*F3(K)-1 371249.0* 
IF 4 (K))) + (3953664. 0*F2 CK) +277020 ,0*F5 (K) ) ) 

RF'lURh 

ELD 



APPENDIX - B 


LIST OF PROGRAM FOR CHAPTER -V 



appendix - V 

************ 


MAIN PROGRAM 

FOURTH PORDER FINITE DIFFERENCE TECHNIQUE. 


Y"(X) + (P1 (X))*Y'tX)-Q(X)*Y(X}=F{X) 
Y'(0)=0 ; YC1 )=BETA. 

WHERE P1(X)=K/X . 


arguments 

********** 

N1 - NUMBER OF MFSH POINTS. 

H - STEP SIZE 

PI, 0 AMO F — THE COFF .FUNCTIONS FROM THE GIVEN 
equation. 

DPI l AMD DDP11 — THE FIRST AND SECOND DERIVATIVES OF Pi, 
DO AND DDQ — THE FIRST AND SECOND DERI AVATIVES OF Q 
DR AMD DDR — THE FIRST AND SECOND DERI AVATIVES OF R 
Ml — MESH PONT WHERE THE SOLUTION TO BE OBTAINED. 


HERE THE INITIAL CONDITION TO BE GIVEN 
BETA=5.5 


DIMENSION XI (0; 100) ,X2C0: 100) ,X3C0: 100) ,Y1 CO: 100) ,Y2CQ:100) , 
1 Y3C 0:100) ,7.1(0;100),Z2(0:100),Z3CO:100),LCO:100), 

1 W(0:100),T(0:100),M(0:100) # N(0:100), 

l PC0:l00),u(0:l0O),vc0:l00),YC0;lO0) 

REAL M , N , L 

********************************************************** 


OPEN CUNIT= 10 , DEVICE* 'USK', PILE* *FD. DAT'} 
R F A D ( 1 0 , * } N 1 , K , M 1 
H=l./FLOAT(Nn 
Y 1 (U)=R( 0 )/( 1 +K) 

Y2(0)s=DR(0)/(l+K) 

Y3(0)*DDR(0)/(1+K) 

Zl(0)=F(0)/(l+K) 

Z2C0)sDF(O)/C1+K) 

Z3CO) = DDFCO)/( l+K) 


M . (0)=2.0+( (H**?)/6)*(b.o*Yl (0) ) + ( (H**4)/12)*Y3C0) 
Tj(0) = t . 0 * C (H**2)/12)*Y1 (0)-((R**3)/24)*(2.G*Y2(Q)) 
f HO) = l ,0-C C (H**2)/!2)*YI (0 ) ) + ( (H**3 ) /24 ) *(2 .O*Y2(0) ) 
P(0)s(H**2)*Zl(0) + ( (!!**4)/12)*Z3C0) 


nu 111 1=1, HI 

m=x+i*u 

X1(1) = P1 ( AI ) 

X2(I)=DP1 (AI) 

X3 ( I ) “PiiPl ( AT ) 

mn=p(Ai) 

Y2(I)=DR(Al) 

Y3(I)=DPH(AI) 

Z1 (I)=F( AI) 

Z2(I)*DFCAI) 

Z3(I)=DDF(AI) 

CONTINUE 

00 20 1=1 ,M-1 

M(I)=2.0+((H**2)/6)*(5.0*Y1 a) + Xl(I)*Xl(I)+X2(I)) + 

1 C(H**4)/12)*CY3(I)+X1 (I)»Y2CT)-X2CI)*Yia)) 



s»a° 

LCI)=1 .0+(H/2)*Xl(I)+((H**2)/12)*(Xl(I)*Xl CI)+X2CI)-Y1 <13 )+ 

1 ( (H**3 ) /24) * f X3 (I )-2.0*Y2 { I) -XI (I)*Y1(I)) 

N(I)*1.0-(H/?)*X1 (T) + CCH**2)/12)*<Xia)*Xl(I)+X2(I)-YlCl))- 
1 ( (H**3)/24)*(X3(I}-2.0*Y2(I)->Ci CI)*Y1(I)) 

P(I) = CH**2)*Zim + C (H**4)/12)*fZ3(I) + XHT)*Z2(I)-X2(I)*ZlCI)) 

continue 


THE INITIAL CONDITIONS FOR W'S AND T'S 
W(0)=CL(0)+NC0) )/M(0) 

T(0)=-P(U)/MC0) 


PU 30 I=i ,N1 -1 

w ( i ) = a . ( i ) ) / f m c i ) - w 1 1 j * w ci' - 1 ) ) 

T(I) = CP(I)-Nm*TCI-l))/(NCl)*W(I-J 

CONTINUE 


Y'CNDsBfcTA 

DO 40 I.=N 1-1 #0.-1 

Y(I)swCn*Y(I + l)+TCI) 

CONTINUE 


TYPE 900 

FORMAT ( 5 X , ' — ' t ) 

DO 45 IsO.Nl ,M1 
TYPE 901 , Y ( I ) 

F0RMAT(5X,E20.10/) 

CONTINUE 
TYPE 902 

FORMATC5X, ' — — — — — — ' 

GO TO 211 



STOP 

END 


SAMPLE PROGRAM 
PROBLEM 5.1 


REAL FUNCTION PlOO 

P1-2.0/X 

RETURN 

ELD 

REAL FUtiCTIOT p(X) 

R=4 .0 

RETURN 

EH.) 

REAL FUNCTION F(X) 
F =«2 • 0 

RETURN 

END 

REAL FUNCTION DPI (X) 

DPt=-2,U/X**2 

RETURN 

END 

REAL FUNCTION DR(X) 

OR=0 . 0 
RETURN 
END 

REAL FUNCTION DFf X) 

r»F*0.0 

RETURN 



EM) 

REAL FUNCTION DDP1 (X) 
r)DPt = 4.0/(X**3) 

RETURN 

END 

REAL FUNCTION DDR (X) 

DDR=0.0 

RETURN 

END 

REAL FUNCTION DDF ( X ) 

DDFsO.C 

RETURN 

END 



APPENDIX - C 


LIST OF PROGRAM FOR CHAPTER - VIII 


APPENDIX - VIII 
*************** 


HAIM PROGRAM 

CUBIC SPLINE METHOD FOR 

INFINITE INTERVAL PROBLEMS • 

*************************************************** 
Y"CX) + P(X)*Y'(X)+Q(X)*Y(X) = R(X) 

YCajaalpha ; Y(INF)=BETA 
WHERE INF DENOTES THE INFINITY. 
*************************************************** 
P, V AND R - THE FUNCTIONS FROM THE 
GIVEN EQUATION 
5! - STEP SIZE 

XT. HE -n - THE FINITE POINT 

APZHRO -ALPHA ZERO , BE ZERO - BETA ZERO, 

AND «\ APE THE CONSTANTS EVALUATED FROM 
'FTP. ASYMPTOTIC BOUNDARY CONDITIONS 
HI - THE MESH POINT WHERE THE SOLUTION 
TO HE OBTAINED. 


DOUBLE PRECISION A(O:200O) , B CO : 2000 ) , C ( 0 : 2000) ,CKO:2GOO) 
I E to: 2000) , D (0 ! 2000 ),T CO! 2000) ,F (OS 2000), 

1 K (0:2000) ,Y CO:2000),XI COS 2000) , 

1 X2 CO : 2000) ,X3 (0:2000), ALPHt ,BETA1,GAMMA1, 

1 APZERU , BE ZERO , S2 ,S3 , SS , SSI , TT , Til , TT2 , 

1 UU,UU1 ,VV, VV1,YY,YY1 


OPEN ( UHIT 5 15 , DEVICE 5 'DSK * ,FILE- *CU8 , DAT * ) 
HEAD(15,*)M,H,APZFR0,BEZERU,K,m1 



UO 21 1 = 0 , N 
AI=X+I*H 
XI ( I ) =P ( A.I ) 

X 2 ( I ) = 0 ( AI ) 

X 3 CI)=R(AI) 

CONTINUE 
DO 2u 1=1, N 

AC j )=(1 .0-H*O.333333333*Xi(I-l)+H*0,3333333*XlCI)-0.08333333* 
1 XI CI-1 )*XICI)*H*H) 

Continue 

or* 3 u i=i, o-i 

bm = U ,0-b*o.3333333*Xl (I )+H*0.3333333*Xl(I + U-Q. 08333333* 
l XI (3.+ n*Xl (I)*H*H) 

CU n = Ct ,( ;> + 0.291bb66*H*(XiCI + t )-XtCT-l) ) -0.08333 33*rt*H* 

i xm-i)*xm+n) 

Tm:*,A(l J, BCI), Cl CD 

CONTINUE 


DO 40 1=1, 0-1 

ea) = (l .0+n*0.5*Xl(I + l)+0.16fj6b66e>*H*H*(X2<I + l)) )*Ht) 
IHI) = U.0-H¥0.5*Xl(I-l)-K>.lt>66ob6b*!T*H*X2CI-l))*BCi:> 
Ca) = Ci ,0+h*0,5*Xi Cl+l))*A(t} + (i.0-ri*Q.5*XHI-l) 

1 C0.66666o6*H*H*X2CI)*CUI)) 

KlI) = -i*H*C .Jbbtobbt»&*CA(I>*X3(I + m4.0*Cl(I)*X3(I) + 

J RCI)*X3(I-.l)) 

TYPK* # Dcn,cCi),mi) 

CUfiTIfUJK 


TO FIND THE INITIAL CONDITIONS FOR W*S AND T # S 



SS=C. 3333 333*X 1CN) -0.166666667*X1 (N)*X1 (N-l 1+0,055555556* 

1 H*H*X1CM)*X2CM-1) 

SS 1-CSS* APZ£R0)/A(N) 

TTs0.083 3 33 333*H*X1(N)*X1(N-1 1+0.1 6666667*X1(N-1 1-0,16666667 
1 *H*X2CN-l)-0,055555556*ri*H*Xl(Nl*X2CN-l ) 

TTls(TT*APZKRO)/A(K) 

TT2=-C (SS1+TT11-APZER0/H1 
TYPE*, SSI ,TT1 ,TT2 

«HI=0. Ill mill*H*H*Xl (N-1)*X2 CM J-0.3 3 333333*X2(Ml-0, 33333 37* 
1 XI 00 + 0,1 6666667 *X1 (N-l )*X1 CO 

ttUl = ( |JU* AP'ZKHO ) / A (N ) 

VV=-i0.1 o6n<>667*XUi4-n+0,08 3 333333 3*H*Xl CN 1 *X 1 (N- 1 )+ 

1 0 , 0 2 7 7 7 7 7 7 8 * H + H * X 1 CN-J) *X2 C fti ) ) 

V V 1 s ( 7 V* APZERO ) / A CH ) 

YY=(UUl+VVl+APZEPn/H+BEZERUl 
■nff»R*,UUl ,W 1 ,YY 

YYl = fC+U ,ltllllll*APZERO*H** 2 *X1CN-1)*X3(N)/B( N-l) -0.3333337* 
1 APZFJF0*H*X3Ctl)/B(N-l l+0.05555556*APZERO*H**2*Xi (N 1 * 

1 X3CH-1)/BCN-1 J-0. 1666667 *APZER0*H*X3CN-1 3/ACN)- 

1 U,05555556*X1 CNl*X3(H-i 1*APZER0**2/ACN 1-0. 027777778* 

l H**2*X1(N-1)*X3CN)/A(N) 

TYPE* , YY1 


AljPHl aYY 
BETA1=TT2 
GAf« •tAlaYYl 


TO STORE THE VALUES OF W'S AND T'S BY USING 
THE INITIAL VALUES OF W'S AND T'S 



W(N-i)=RETAl/ALPHl 
00 103 1=1, N-l 
0 1 =M~ 1 
J=J1-1 

w(j)=ocji)/(C(jn-Ecan+«(ji)) 

CONTINUE 

T(N-l)=GAMMAl/ALPhl 
DO 104 1=1, Si-1 
K 1 = M- I 
K=K1 -1 

TCK) = CF(n )*T(K 1 )-F (Fl))/CC(Kt)-E(Kl)*WCKn) 

1: 1 wue 


• IKRfc t'Hg ISUHAL CONDITION' TO BE GIVEN. 
V ( U J = 1 . 0 - 


USING THE VALUES OF W(t)'S AND LCD'S THE 
SOLUTION yen's ARE calculated by forward 
PROCESS 


or- 90 1 = 0, N-l 

mo ):.sa)oii)i!(i) 

CX?E*,I,Y(1) 

coot r hue 

TYPE 210 

FORMAT C 12 A, ' */) 

I t PC 212 

FORMATC13X, 'COMPUTED SOLUTION'/) 

Di) 1 0 U 1 s0 ( n, #1 
TYPE 213, VCD 



FO H M A T ( 1 0 X , E 2 0 , 1 0 / ) 

c n h riNUF 

TYPE 214 
FORMAT (12X, ' — 

GO TO 111 

STOP 

END 


SAMPLE PROGRAM 
PROBLEM 8.1 


REAL FUNCTION PtX). 

P=2 * C 

RETURN 

END 

REAL FUNCTION OtX) 

y=-2.o 

RETURN 

END 

REAL FUNCTION PCX) 
R=-F.:XP(-2,0*X) 



fY) fr-T H - t S> - 0- Ai- A/t'/V) 



